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FOREWORD 


This two part final report summarizes the efforts and accomplishments of 
the second phase of the subject contract covering the period June 1975 through 
September 1977. The study was conducted for the NASA-Lewis Research 
Center, Cleveland, Ohio, by personnel in the Engineering Sciences Section, 
Lockheed -Huntsville Research & Engineering Center, Huntsville, Alabama. 

The NASA-LeRC Project Engineer was Dr. C. C. Chamis, Mail Stop 49-3. 

C. H. Lee was the principal investigator for the study, and the study was 
supervised by B.H. Shirley. 

Work was begun on this contract on 28 June 1974. Phase I of the study 
was completed in June 1975 and an interim report was submitted and published 
as NASA CR-134933. Efforts on the study were directed toward the local analy- 
sis of three-dimensional high velocity impact problems based on Eulerian repre- 
sentation. The follow-on study in Phase IX was directed toward the global finite 
element analysis of the problem by coupling the Eulerian mode, based on the 
results obtained in Phase I, with the Lagrangian mode. An interfacing pro- 
cedure for coupling the presently developed program, CELFE, with NASTRAN 
was also accomplished to increase the range of applicability for the program. 
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ABSTRACT 


The theoretical formulation of the Coupled Eulerian-Lagrangian Finite 
Element (CELFE) program for high velocity impact analysis is described. 

The formulation comprises the development of general objective of the pre- 
sent study, a 3-D finite element program capable of simulating the dynamic 
behavior in the vicinity of the impact point, together with predicting the dy- 
namic response in the remaining part of the structural component subjected 
to high velocity impact. 

The finite element algorithm is formulated in a general moving co- 
ordinate system. In the vicinity of the impact point contained by a moving 
failure front, the relative velocity of the coordinate system will approach 
the material particle velocity. The dynamic behavior inside the region will 
thus be described by Eulerian formulation based on a hydroelasto-viscoplastic 
model. The failure front which can be regarded as the boundary of the im- 
pact zone is described by a transition layer. The layer will change the repre- 
sentation from the Eulerian mode to the Lagrangian mode outside' the failure 
front by varying the relative velocity of the coordinate system to zero. The 
dynamic response in the remaining part of the structure described by the 
Lagrangian formulation can thus be treated using advanced structural 
analysis. 

An interfacing algorithm for coupling CELFE with NASTRAN is con- 
structed to provide computational capabilities for large structures. The 
coupled CE LF E/NASTRAN system also covers the computations of the struc- 
tural response for a large time after the impact process is completed. Locally, 
the system has provisions for projectile rebound and for projectile sliding due 
to oblique impact. 
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1. INTRODUCTION 


Over the last two decades considerable interest has been shown in the 
study of high velocity impact problems. This is primarily because of the 
extensive development of space vehicles and the need for information con- 
cerning the meteoroid hazard to these vehicles. Current interest in high 
velocity impact studies is largely due to the concern over the impact of inter- 
planetary debris on space vehicles. These research results also have other 
industrial, as well as military, applications; for instance, in the high velocity 
impact on turbine blades of foreign objects such as rocks, birds or metal 
projectiles. 

Motivated by military applications, seven symposiums were held on 
the subject during 1955-65 [1-7] under the sponsorship of the Army, Air Force 
and Navy. The bulk of the material presented at these symposiums consisted 
of experimental results and very little appeared on the theories explaining the 
high velocity impact phenom ena. During the late 1960s efforts were expended • 
toward formulating realistic theories to explain the complex dynamic and 
mechanical response of materials in hypervelocity impact [8-12] . 

The term "hyper velocity 1 ' (or high velocity) refers to the impact velo- 
city regime in which the maximum stress developed by the primary shock 
wave greatly exceeds the material strengths of both the target and the pro- 
jectile. During the early stages of the impact process, the target and the 
projectile behave essentially as compressible fluids and, consequently, 
several researchers [13-16] have employed pure hydrodynamic models to 
analyze the hypervelocity impact problems. However, these shock stresses 
decay very rapidly as the wave propagates away from the impact point and 
reach values comparable to the material yield strengths. From this point 
onward the material strengths become important in determining the stress 
and velocity fields in the target material and, consequently, due to more 
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involved material response to these rapidly applied stresses, the process 
gets more complicated. Hence, the purely hydrodynamic model cannot ade- 
quately describe the complex physical phenomena [17] and a more realistic 
model must be employed to account for all aspects of the high velocity impact 
phenomena. Several researchers have included the strength effects (e.g., 
see [18]) and there appeared many computer codes using more realistic 
hydrodynamic-elastic-viscoplastic models, [10] and [19-23]. The numerical 
techniques were largely based on the finite difference method, and they appear 
to be suited for problems- with relatively simple geometries. 

Since a hydrodynamic -elastic -viscoplastic model is basically a structural 
model, it would be appropriate to use structural techniques to solve the problem. 
Also, during the impact process, the geometry involved is generally very com- 
plicated, which would necessarily require a versatile and flexible technique in 
order to treat it accurately and realistically. The finite element method, which 
has proved to be highly successful in analyzing structural problems, is considered 
such a candidate. Leimbach and Prozan [24] have shown the superiority of the 
finite element method over the finite difference method for a simple impact 
problem, although their model is rather simple and unable to predict the actual 
dynamical response of materials in the impact region. 

The general objective of the present study is to develop a three-dimensional 
finite element code to analyze structural components subjected to high velocity 
impact, using governing equations based on Eulerian and Lagrangian formulations. 
In particular, the impact point with its vicinity is td be represented by an Eulerian 
hydrodynamic -elasto-viscoplastic model, while the remaining structural com- 
ponents are to be analyzed with existing computer programs such as NASTRAN, 
for structural analysis. To bridge the gap between the Eulerian mode used for 
the impact region and the Lagrangian mode generally employed in structural 
analysis, a finite element code based on coupled Eulerian and Lagrangian form- 
ulation is developed. 
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This report contains the development of the complete system of the 
Coupled Eulerian-Lagrangian Finite Element (CELFE) program, including 
the interfacing procedure with NASTRAN, for handling a given structure sub- 
jected to high velocity impact. The report is divided into two parts: Part I 
includes a mathematical model for the impact process, together with the 
theoretical aspects of the CELFE algorithm. Part II presents listings and 
user's manual of the CELFE program and an illustration of the CELFE/ 
NASTRAN interfacing procedure. 

Following this introduction, Section 2 will present a preliminary dis- 
cussion on the mathematical model to be developed, the basic idea of the 
associated simulation procedure, together with the geometric description of 
the Eulerian, Lagxangian and the coupled Eulerian-Lagrangian- formulations . 
According to the simulation procedure, the geometry of the entire structure 
subjected to high velocity impact is divided into two parts, namely, the im- 
pact zone covering the impact point and its vicinity, and the Lagrangian 
(structure) zone for the remaining part. In the impact zone, the dynamic 
behavior induced by the impact is described by 'a hydroelasto-viscoplastic 
model represented in a general moving coordinate system. The development 
of the model and the associated governing equations are discussed in Section 3. 
By assuming the small deformation theory, on the other hand, the dynamic re- 
sponse due to the impact in the Lagrangian formulation is predicted using 
classical structural dynamic analyses. Sejtion 4 summarizes briefly the re- 
lated elastodynamic equations. In Section 5, a finite element algorithm for a 
general system of quasilinear equations with conservation form is formulated. 
The finite element analog of the equations is constructed as a direct conse- 
quence of the theorem of weak solutions. Thus, the entropy condition can be 
satisfied automatically in the formulation . The algorithm is formulated 
mainly for simulating the viscoplastic flow in the impact zone. In the same 
section, a generalized two-step scheme for time integration is also proposed 
to remedy the numerical instabilities and at the same time to speed up the 
computations. This scheme is applied for the entire CELFE procedure. A 
description of the complete CELFE procedure, the mechanisms coupling the 
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Eulerian- Lag rangian modes in particular, is presented in Section 6. The 
coupling mechanisms consist of an in-core coupling procedure to smoothly 
translate the coordinate system from Eulerian into Lagrangian, and an inter- 
facing procedure for the CELFE/NASTRAN system. The complete CELFE 
system covers the computations for large time after the impact process is 
completed; it also simulates the rebound of the projectile from the target, 
ahd the sliding due to. oblique impact. The criteria for the various cases are 
also discussed in Section 6. Section 7 summarizes results from numerical 
experiments for various test problems. The test problems include an inviscid 
hydrodynamic model with Eulerian description; a hydroelasto -vis coplastic 
model in Eulerian and the general moving coordinate system; and a complete 
CELFE model. Section 8 summarizes the present study and describes pos- 
sible improvements and extensions of CELFE as related to computational 
efficiency, capability, and application. 
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2. A HYDROELASTO -VISCOELASTIC MODEL 
FOR HIGH VELOCITY IMPACT ANALYSIS 


The high velocity impact of a projectile with a solid target results in 
extremely complex phenomena. A complete description of this problem would 
involve considerations of all phases in the theory of continuum mechanics. 

This includes not only the compressible fluid flow, dynamics of elasticity and 
plasticity, but also other behavior such as melting and solidification, vapori- 
zation and condensation, and the kinetics of phase change. As a consequence,, 
certain simplification is needed for tackling the problems. Several models 
have been proposed for various stages of the problem (see, e.g., [13, 16 ], [19- 
23]). In this section, a brief introduction of a hydrodynamic-elastic-viscoplastic 
model coupling the Eulerian and Lagrangian modes is presented. 

2.1 GENERAL DISCUSSIONS OF THE IMPACT PROCESS 

The analysis of high velocity projectile mechanics can be divided into 
two parts: (1) dynamic behavior of the projectile and target during the pene- 
tration process, and (2) structural response of the target after the' penetration 
process is completed. The complex phenomena occur mostly during the pene- 
tration process. The following discussions are directed primarily toward 
this stage. 

During the short period of time in which the projectile contacts the 
target, a plane shock wave is generated in the projectile as well as in the 
target. The pressure and hence the temperature behind these shock fronts 
are the highest that exist throughout the entire impact process, and may 
cause material failure in both the target and the projectile in the vicinity 
of the impact point. 

Immediately after the impact, due to the high pressure, the material 
strength of the target can be ignored, and the material can be assumed to 
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behave essentially as an inviscid, compressible fluid. The shock waves gen- 
erated in the projectile and target travel away from the interface (see Fig. 2-1). 
If the projectile is finite (in diameter), rarefaction waves will be generated and 
transmitted toward the axis of symmetry. 

Formulation of rarefaction waves, results in the ejection of target and 
projectile material particles. Moreover, the rarefaction waves weaken the 
shock wave and change its plane shape to approximately spherical. The 
strength of the shock continues to decrease due to the spherical attenuation 
and additional rarefaction waves, and the influence of the material strength 
and the strain rate effects must be taken into account. At the same time, the 
failure region grows following the outgoing shock until the penetration process 
is completed. In this stage, the dynamic behavior in the impact zone covering 
the failure region can be considered practically as a compressible, non- 
Newtonian flow; and the structural response in the remaining part of the target 
begins to appear as the disturbances propagate through the material. 

2.2 HYDROELAS TO -VISCOPLASTIC MODEL 

Based on the discussions presented in the previous subsection, an ani- 
sotropic, hydroelasto-viscoplastic model for high velocity impact analysis 
can be set up as follows; 

During a certain time period covering the penetration process, the dy- 
namic behavior in the vicinity of the impact point containing a failure zone 
with moving front is regarded to be characterized by a compressible flow. 

The flow is governed by the Navier-Stokes equations (i.e., the conservation 
equations of mass, momentum and energy), in conjunction with an appropriate 
equation of state and the constitutive equations for both the target and pro- 
jectile materials. In addition, a system of well established [25] , hydrodynamic - 
elastic -viscoplastic constitutive equations based on Prandtl-Reuss theory is 
adopted to describe the material properties. 
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Fig. 2-1 - High Velocity Impact Process 
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If there is no substantial heat flux across the boundary of the body, the 
possible heat transfer will be induced only by heat source generated during 
the. impact process-. This heat transfer is - generally not important because; 

m There is a large temperature gradient across the failure front; 

9 A large part of the heat produced during the impact process, inside 
the failure zone, is converted to latent heat absorbed by the ma- 
terial for possible phase change. 

Thus, at the present stage, the heat transfer will not be considered. When- 
ever the thermal effects between the body and its environment become impor- 
tant, these effects may be added into the energy equation in a straightforward 
manner. 

The phase change, on the other hand. Is assumed to be included into the 
equation of state as an approximation. There will be no attempt at this stage 
to incorporate the diffusion equation governing the phase change into the im- 
pact problem due to the above mentioned complexity. From a practical view- 
point, an accurately correlated, empirical equation of state would provide a 
better result on accounting the phase change in the present stage of develop- 
ment. 


For three-dimensional problems, the formulation consists of five con- 
servative equations and nine constitutive equations that are nonlinear partial 
differential equations in space and time. These equations will be solved using 
the finite element method in conjunction with an appropriate equation of state. 
The failure of the material caused by the impact is computed by various cri- 
teria depending on the material as well as the numerical scheme. In the pre- 
sent task where the failure zone is simulated using the Eulerian mode, criteria 
for the onset of viscoplastic flow are used. For composite materials, the 
failure prediction follows Chamis' criterion developed using micromechanics 
structure theory [26-28] ; for metals, the plastic yielding is predicted in accord 
ance with von Mises criterion. 
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Away from the impact zone, however, the structural response of the 
material is expected to be small compared to the one in the impact 'zone, and 
can be described essentially by the classical theories in structural dynamics. 

If the above mentioned system of 14 equations were to be retained for com- 
puting the dynamic response in this region, some of the equations become 
practically redundant, and thus a massive computation. which occupies a large 
computer storage turns out to be trivial. In view of this inefficiency^ the 
classical elastodynamic analyses are adopted for this portion. 

The governing equations in this part depends on the choice of the ana- 
lytic techniques. They may vary from the displacement, stress (or force), and 
energy finite element method, to the mixed variables and hybrid finite ele- 
ment method. The choice of analytic technique follows not only the type of 
structure under consideration, hut it is also required to meet the following 
objectives: it must have the computational efficiency and simplicity to couple 
the analysis with that in the impact zone (i.e., in coupling the Eulerian and 
Lagrangian modes); and it must be capable of incorporation with available 
structural analysis programs such as NASTRAN. 

.2.3 COUPLED EULERIAN -LAGRANGIAN REPRESENTATION 
2.3.1 Geometric Considerations 

Consider an open (i.e., not including its boundary) bounded region 
at a time t = 0 in three-dimensional Euclidean space with its boundary 9<T2 
(see Fig. 2-2). The union of and its boundary is the complete region 
(occupied by a' body at time t = 0) and is denoted by 

9 Q = n Q U 80 

where U denotes the union of two sets. In , there is a region £2^ with moving 
boundary Let denote the complement of 9 ^ inf2 Q . Note that the region 

^2 has two boundaries: one external boundary 3f2 and one internal (the interface) 
boundary 802^. 
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Fig. 2-2 - A Three-Dimensional Region at t = 0 

Suppose that we are given, at time t = 0, the state of the fluid (or solid) 
occupying the region the external forces acting on the boundary, and the 
boundary conditions. We are then required to determine the state of the fluid 
and the shape of the region at a subsequent time, t = T. Since we plan to 
analyze the impact point and the surrounding region, £1 separated from the 
remaining region, it is convenient to consider two different fluids (or ma- 
terials) occupying two different portions, £2 ^ and of that is 

The surface is the material interface of the two regions, and it moves, as 
time advances, in such a manner that there are jumps on the pressure (or 
stresses) and velocity components over the surface. 

According to this specification, we may define: 

• Q j to be the Impact Zone which contains the impact point and the 
failure zone with a moving boundary 9 q and 

• i° be the Structure Zone. 
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2.3.2 Euler and Lagrange Descriptions 


There are two basic descriptions well known in continuum mechanics 
with respect to which the governing equations can be derived and computa- 
tions can be carried out. In the Lagrange description, mostly used in solid 
and structural problems, a coordinate system is fixed in the body or con- 
figuration £2 to be studied. The deformation of the projectile -tar get con- 
figuration is then measured with respect to this coordinate system. 
Consequently, the positions of the boundary 3J2 and of the material interface 
are automatically determined. The description also permits the use of 
constitutive relations for the material in which the stress history of each 
portion of the body is taken into account. However, the Lagrange description 
is totally unsatisfactory for calculating a large deformation, and a flow in 
which turbulence develops or in which new material interfaces develop. In 
this case the nodal points of the mesh will attempt to follow the motion and 
the material particles which were initially adjacent to each other in £2^ no 
longer remain so. 

In the Eulerian description, mostly used in fluid mechanic problems, the 
coordinate system is fixed i n the space rather than in the body or configura- 
tion, and the calculations follow the material point that happens to.be in a 
given location at that particular time. In this case the large distortions do 
not cause any problems; however, more than one material cannot be treated 
accurately. The curves approximating 3$2 and 3f2^ move with the body and 
therefore create irregular, time -dependent boundary zones in the fixed 
Eulerian mesh. 

2.3.2 Coupled Eulerian-Lagrangian Description 

Obviously, neither Lagrange description nor Eulerian description alone 
is ideally suited for the analysis of impact problems. It is both natural and 
desirable to combine the use of Eulerian and Lagrangian modes depending on 
whether the material is in a fluid state (Eulerian mode) or solid state (La- 
grangian mode). In doing so, there will be a great deal of flexibility in ap- 
proximating the problem, thus enhancing the solution process regarding 
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accuracy and computational efficiency. The idea is very similar to the well 
known sub structuring technique, but now with different descriptions (or modes) 
used' in various regions. More specifically, in analyzing the impact problem 
os structural components, the existing NASTRAN program (in Lagrandian 
mode) can be used to advantage for the structure part, while an Eulerian 
description is necessary for the impact point and its vicinity. 

The basic idea of coupling fhe Eulerian and Lagrandian modes is to 
approximate the structure in consideration, £l Q > into the combination of 
Eulerian and Lagrangian subregions separated by a transition zone. Figure 
2-3 shows a configuration consisting of these different subregions. 

As shown in Fig. 2-3, the structure zone ,Q.,, can be identified with the 
Lagrangian subregion. Thus, it may also be called the Lagrangian zone , and 
denoted as L zone . In this subregion, the problem will be formulated accord- 
ing to the classical theory of e las tody namics under the Lagrangian coordinate 
system. 

The impact zone , oh the other hand, consists of the Eulerian subregion 
and the transition zone. We denote the Eulerian zone as E zone, and the tran- 
sition zone as E-L zone represented by a moving coordinate system having 
varied local velocity with respect to the material particles. The relative 
velocity varies in such a way that, for each fixed t, it becomes identical to 
the local velocity within the E zone, and decreases monotonically toward the 
moving boundary ^ through the E-L zone, and vanishes on This im- 

plies that the moving boundary ^ is now described by the Lagrangian system. 

The formulation thus enables us to trace the propagation of the failure 
front at each t on one hand, and couple the Eulerian mode in the E zone and 
the Lagrangian mode in the L zone smoothly for all t > 0 on the other hand. 
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Fig, 2-3 - A Typical Configuration, in Coupled Eulerian-Lagrangian 
R epx e s entati on 
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3. BASIC EQUATIONS AND PROBLEM FORMULATION 
IN THE IMPACT ZONE 


In this and the following section, the ha sic equations for various zones 
based on the model discussed in Section 2 will be formulated. As mentioned 
previously, these equations may be composed of different sets of primary 
variables. A prerequisite of choosing the equations for each zone is that 
there must exist at least one set of primary variables, i.e., the coupling 
variables , either common or directly related to the ones in the remaining 
zones. Since the final objective of this study is to couple the CELFE pro- 
gram with the NASTRAN program in which the displacement vector is used 
in the dynamic analysis of structures, the displacement will be chosen ex- 
clusively as our coupling variables in developing the CELFE code. 

To couple the Eulerian and Lagrangian modes smoothly as outlined in 
Section 2, it is more convenient to represent the governing equations for the 
entire impact zone in a general moving coordinate system, instead of using 
the Eulerian mode in the E zone and the moving coordinates in the transition 
(E-L) zone. Mathematically, this representation can provide an identical 
approximation on changing the system from Eulerian to Lagrangian for all 
time as the finite element approximation of the physical problem itself. 
From a computational viewpoint, this representation cannot only simplify 
the programming, but can at the same time handle the movement of the free 
surface automatically. 

3.1 LOCAL REPRESENTATION OF DYNAMICAL SYSTEMS UNDER 
AN ARBITRARY FRAME 

Let x = (xj.x^x^) be a Cartesian coordinate system moving with a 
certain velocity relative to the motion of the physical field, and let X = (X^, 
X^, ^ 3 ) k e the reference system. Then, the two coordinate systems can be 
related as 
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x = x(X,t) 


(3.1) 


When we know such a relation for every material particle in the medium, we 
say that the history of the motion is known at time t. We assume that the 
motion is continuous, single valued and that (3.1) can be inverted to give the 
initial position or material coordinates, X. A necessary and sufficient condi- 
tion for the inverse to exist is that the Jacobian should not vanish: 


J = 


3x. 

i 

ax. 

j 


/ 0, 0 < J < oo V" t > 0 


(3.2) 


Let Q denote the velocity of a material element relative to the x 
system, i.e., 


d x 

£2 as = V - X 

^ dt ^ ~~ 


(3.3) 


where is tjie velocity of the material element, and x = 8x/8t is the local time 
rate of change of the x system. If F represents a certain physical quantity 
satisfying a dynamical syst em , then the total time derivative follows the equality 


dF _ 3F 

dt = ‘ at 


+ (n • v) f 


(3.4) 


Clearly, when x = 0, the systems x and X are identical, and x represents the 
Eulerian system: when x = v , on the other hand, it represents the Lagrangian 
system. 

Following the implicit differentiation rule, the divergence of the local 
rate of change of x system can be related to the Jacobian as 


V «x = j/J 


(3.5) 
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thus, (3.3) yields 


V«Q = V - v - j/J 

rV> ' 


(3.6) 


Suppose now that there is a dynamic system with physical quantity F 
satisfying the equation 

= G(x,t) (3.7) 

In particular, the continuity equation in x system is readily obtained as 


+ PV-Z = 0 '• < 3 * 8 > 

where p is the density of the material. Then, the dynamic system (3.7) can 
be deduced by applying (3.3), (3.6) and (3.8) into the following form: 

P.teZi sfrpFA + V .(gp F) = p G (3.9) 

3 t 

where 

A = j/J (3.10) 

It is clear that the dynamic system (3..9) would reduce to the Eulerian form 
when x = 0; and to the Lagrangian form when £2=0. 

3.2 CONSERVATION EQUATIONS 

Let F, defined In Eq. (3.9), be 1, v, and e 4 then the conservation of 
mass, momentum and energy can be obtained and written, respectively,. in 
tensor form as 
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6P | 8<n fc P> 

9t 9x, 


+ p A = 0 


av. 8(« k v-) 

1 A J. + y 

r\j_ r r\ 1 v * 


Bov 


at 9x 


A = 


& 


a + E. 

) 


as + !f^ +EA - 

at ax, ax, 


+ v. F. 
3 3 


for j, k = 1, 2, 3, where 


(3.11) 


(3.12) 


(3.13) 


V. = pv., E = pe 

denote the momentum and total energy, respectively, o\j, i, j = 1, 2, 3 
denote the stresses, and F- now denotes applied or body force. 


When fi. = v., j = 1,2,3, A=Q. Thus, Eqs . (3.1 1) through (3.13) reduce 
to the well-known continuity, momentum and. energy equations, respectively, in the, 
Eulerian description with conservation form: 


9t 


+ 


3( y k P) 



(3.14) 


9V. 9 ( v k V ^ aor., 

1 + iU_ = ___jk +F 

9t 9x^ ' 9x k j 


(3.15) 


9E 

at 


' 9(v k E) 


9x, 


, S(Vj P jk ) 


9x, 


+ v. F . 

1 J 


(3.16) 


When . = 0, 
representation: 


j - 1, 2, 3, Eqs. (3.11) through (3.13) become the Lagrangian 
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f-+pA = 0 


(3.17) 



+ V. A 
3 



+ F. 

3 


9E 

at 


+ E A 


9(v, cp ) 

— ^ — + v. F. 
3x v 3 3 


Here, Adefined in Eq.(3.10) yields 

j - p c /p 


(3.18) 

(3.19) 

(3.20) 


which, is clearly equivalent to the continuity equation (3.17). 

3.3 CONSTITUTIVE EQUATIONS 

When the medium under consideration is inviscid, there are no shearing 
stresses. The stress tensor then becomes 


— cr.-. = -Pi. 
. 13 ij 


where P is the the thermodynamic pressure, and 8„ is the Kronecker delta. 
In general, the stress tensor cr. . can be decomposed into a (dynamic) pressure 

13 

term and viscous terms that depend on the stress versus strain-rate relation. 
Accordingly, the stress tensor and the associated strain rate tensor are re- 
lated to their respective deviatoric tensors by the following forms: 


S.. 

13 


cr. . 
J-3 


- 1/3 a 


kk 


8 .. 

13 


(3.21) 



d tj - ^ “kk 6 


13 


( 3 . 22 ) 


where S.. is the deviatoric stress tensor, d. . is the strain rate tensor defined 

13 13 


by 


d.. 

13 


1 

2 


9v, \ 

8x. + 3x. ) 
J 1 


(3.23) 
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and dL is the deviatoric tensor of strain rate. The normal components of 

are directly related to the pressure induced by the impact; the viscous 

effect, on the other hand, stems from-a constitutive relation of S. . toward d.'.. 
’ ij n 




3.3,1 Deviatoric Stresses -Pres sure Relation 

Let us denote the mean of the principal components of stress by 

Q e - 1/3 c kk (3.24) 

The quantity Q is called the dynamic pressure. The rate of dynamic pres- 
sure dQ/dt may, in general, be decomposed into a thermodynamical reversible 
component dP/dt and a dissipative component d 7 r/dt, i.e., 


dQ _ dP , djr 
dt dt dt 


(3.25) 


In general, the rate at which the local thermodynamic equilibrium is attained 
is much greater than the rate at which a disturbance can be propagated. It is 
then reasonably accurate to assume that the local thermodynamic equilibrium 
exists at each instant. Hence, the reversible rate of pressure, dP/dt, is not 
path -dependent, and its integral P follows the equation of state which can be 
expressed as a function of density p and internal energy e, i.e., P = P (p, e). 
The dissipative component d?r/dt, on the other hand, is generally a path- 
dependent function related to the bulk viscosity and the volume rate of change. 
It characterizes the physical dissipative rate of dilatation. For isotropic 
materials. 


7T = - (A + 2/3 p) d kk 


(3.26) 


The coefficient of bulk viscosity is defined by 

P' = (A +2/3 p) = (P - Q)/d kk (3.27) 


3-6 



If the volumetric changes of the materials are elastic, the change of dynamic 
pressure is reversible which implies that dQ/dt = dP/dt. Then the path- 
independent nature of P yields 

Q = P = P (p, e) (3.28) 

Noting that for incompressible materials d, =0, Eq. (3.28) follows 
immediately from Eq. (3.27). Equation (3.28) is also true for compressible 
fluids for which 

(A + 2/3 11) = 0 (3.29) 

Condition (3.29), known as Stoke's hypothesis, is a reasonable assumption 
for flow of monatomic gases; however, it is not valid for polyatomic gases 
or liquids, and distinction must be made between the mean stress Q and the 
thermodynamic pressure P. 


For high velocity impact problems, the thermodynamic pressure is 
very high, and as aforementioned, the dynamic response in the material can 
thus be considered as an isentropic process. This implies that the material 
under high velocity impact ca n be assumed to possess elastic changes in 
volume, and (3.28) follows. Thus, (3.21) can be written as 


S.. = <j . . + P 6.. (3.30) 

il i] i] v ' 


3.3.2 Hydrodynamic-Elastic-Viscoplastic Constitutive Relation 


For linear elastic materials, the constitutive equation for stress and 
deformation rate is given by 


S.. 

ij 


E ijki d ki 


(3.31) 


where E^.^ is a fourth-order tensor of material parameters. If the stress 
components o\ . are symmetric, E.. ^ is also symmetric, i.e., = 

E jiik’ etc. In general, there are 21 independent elastic constants. For ortho- 
tropic materials, the number of independent constants reduces to 9, and for 
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isotropic materials it reduces to 2. The isotropic constitutive relation is 
given by 

S. . 2 2,1 d. . * Mkk 6. . (3.32) 

where A and fl are the Lame's (or viscous) constants. 


Unfortunately, there does not exist a constitutive relation that describes 
all aspects (elastic, viscous and plastic) of mechanical behavior in a single 
expression, Constitutive relations of plasticity and viscoelasticity are essen- 
tially dynamical in nature. The constitutive relation in classical plasticity 
involves tensors of stress and rate of deformation, and describes rigid, per- 
fectly plastic behavior. When elastic effects are to be considered in the 
analysis, this relation applies to the plastic part of the rate of deformation 
tensor. Similarly, the constitutive equation for viscoelastic material involves 
stress, elastic rate of deformation and the rate of deformation of a viscous 
fluid. In both cases, the elastic component of the rate of deformation tensor 
is usually written as a function of stress rate. Although the stress tensor is 
objective (axiom of objectivity requires that if a stressed continuum performs 
a rigid body motion and the stress field is independent of time when referred 
to a coordinate system attached to and moving with the material, the stress 
rate must vanish identically), the stress -rate tensor is not. Therefore, the 

stress rate do\ ./dt should not occur in this form in the constitutive relation. 

^ J 

An objective tensor containing the stress rate must be defined. Several ob- 
jective tensors containing the stress rate can be constructed. One such tensor 
is due to Jaumann (see [29-] ) and is given by 

a __ 

°k S. ~~ dt + °km U mi “ °mi W km (3.33) 


where 



(3.34) 
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The tensor <7^ is called a stress flux . Other stress. fluxes may be obtained 
by adding objective tensors such as right -hand side of (3.23): 



do. 


k £ 

dt 


+ a .v . 4- a, v . 
ml m, k km m, £ 


do. , 

a _ k£ v 

0 k£ _ dt °m£ k, m- °km i,m 


(3.35) 


a d(y k£ 

°k£ ~ dt” ~ <7 m£ V k, m °km V £, m + a k£ V m,m 


The stress flux cr^ measures the rate of change of the stress component with 
respect to a rectangular Cartesian system that participates in the rotation of 
the material, and cn . = 0 implies that the invariants of stress tensor are sta- 
tionary. 


From Eos. (3.21) and (3.34), it is clear that the stress in Eq. (3.33) can 

A a 

be replaced by the deviatoric stress flux S. .. In general, the stress flux cr. . 

a ij 13 

and thus the deviatoric stress flux S„ are balanced by certain forcing func- 
tions $. Symbolically, we have. 


A 

S. . 
!3 


$ (S. ., d. v.; x., t) 
13 13 3 3 


(3.36) 


Here, the functional $ of S.. ., d.. and v. may be regarded as a generalized 
material potential. The form of the potential depends on the material under 
consideration as well as its failure criterion, and may also include the effects 
of phase changes and related properties. 


Assume that. the materials under consideration are linear and follow 
the Prandtl-Ruess flow rule [25 ] - . In accordance with the rule, the rate of 
deformation tensor (strain -rate tensor) d„ can be split into two parts, i.e., 

d.. = df. + d7? (3.37) 

13 13 13 
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where df. is the elastic component and df? is the viscoplastic component of 
the deformation rate tensor. Here it is assumed that the medium is initially 
unstressed. The plastic part of the strain -rate tensor -is- at any instant pro- 
portional to the deviatoric stresses, i.e.. 


df? = A.. 




ijki S ki 


{3.38) 


where > 0 and is in general a function of both time and spatial coordi- 

nates. The elastic part, on the other hand, is related to the stress flux 

A 

S. . as 

13 


S.. = C... . df. 
13 ljki kl 


(3.39) 


Here the coefficient matrix depends on °f 33 q. (3.31). In this 

study we are concerned mainly with the fiber composite (orthotropic) and the 
isotropic materials. 


For generally orthotropic materials, the local stiffne 
can be expressed as 


ss matrix C 


ijki 


C 11 

C 12 

CO 

i—i 

u 

0 

0 

0 

C 21 

C 22 

C 23 

0 

0 

0 

C 31 

C 32 

C 33 

0 

0 

0 

0 

0 

0 

C 44 

O' 

0 

0 

0 

0 

0 

C 55 

0 

0 

0 

0 

0 

0 

C 66 _ 


(3.40) 


where the subscripts a, (3 represent the location in the array with a indicating 
the effect and (3 indicating the cause. Comparing with subscripts ijki , we note 
that 1 » 11, 2 > 22, 3 «— * 33, 4 23, 5 31, 6 12 as a versus ij and J3 

versus ki; and 

C a{3 = C f 3 a 5 or C ijki = C kiij 
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In terms of engineering constants, the complete stiffness matrix with 
respect to the (x^, x^, Xg)-system with x^, parallel to the fibers is given by 


C H ~ ^ H ~ v i 9) H E„, 


23 32 11’ 22 


3 ! ‘' 13 ' 22 ’ 


C 33 ^ v \2 v 2l ) H E 33’ 


C 1 2 (v 21 + ^23 ^31^ H E ll’ C 13 “ (l/ 31 + ^21 v 32^ E 11 


^23 4 9 ^91 ) E99, 


32 12 31 22' 


(3.41) 


C 44 = G 23’ C 55 = G ^' C 4A = G 


13' 66 


12 


H = (1 - "l2 V 21 - " 2 3 "32 - "31 "13 ' 2 "12 "23 "31> 


-1 


where E represents Young's modulus, G the shear modulus, and v the Poisson 
ratio. The subscripts ij in E,G, and v represent the measure of the strain 
in the x. direction under uniaxial normal stress in the x. direction. 

] 1 

For isotropic materials, the stiffness matrix becomes 


C. . = 2 {1 <$.„ + A 6- . 6, „ 

ljki ^ lk jf ij ki 


(3.42) 


where JLl and X are Lame constants. Thus, Eq. (3.39) can be deduced further 
as 


a d S.. 

S. . = — ^ + S. co . - S . co- 
ij dt 1 m mj mj 1 m 


= 2 (1 d. . 

i] 


(3.43) 


Under the general moving coordinate system discussed in Section 3.1 
(cf., Eq. (3.4)), Eq. (3.39) can be written as 


as.. 
II 

at 


as. . 




k 3x, 


11 - 


= C. .. . d, . - S. co . + S . co- 
ljkf ki lm mj mj im 


(3.44) 
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Similarly, by letting £2^. = for all k = 1, 2, 3, (3.44) can be deduced to the 
Eulerian mode as 


3S. . 9S. . , 

at 1 + V k ~ C ijki d kf " S im W mj + S mj ^im 


(3.45) 


and to the Lagrangian. mode by £2^. = 0 as 


3S. . , 

q, 1 ' ^ = C. .. . df. - S. co . + S . 

3t ljkf ki lm mj mj im 


(3.46) 


3.4 FAILURE CRITERIA 

e 1 

Notice that the elastic part of the deviatoric strain -rate tensor, d„, 
in Eq. (3.39) can be related to d). and its viscoplastic counterpart as shown in 
Eq. (3.37). The viscoplastic part, which is in turn related to the deviatoric 
stresses, depends on the failure criteria of various materials. 

3.4.1 Plastic Yield for Isotropic Materials 

For isotropic materials, the plastic behavior is assumed to take place 
when a certain function of deviatoxdc stresses vanishes. This function is 
called the yield function. The yield function is constructed based on the 
following assumptions: 

1. The yield surface is convex (or smooth). 

2. The plastic component of the deformation rate tensor is normal 
to the yield surface at a smooth point. 

If a yield condition is given by 

F(cr *) = f(cr K) - Y(tf) = 0 (3.47) 

1 J IJ 

A 

with F < 0 denoting the purely elastic region, then the deformation rate tensor 
will be a function of positive values of f. Here Y is a yield stress and K is a 
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history dependent hardening (or softening) parameter. We now introduce the 
notion of ''plastic potential" defined as 


= V (f / f o ):> ^- (3.48) 

•ij 

where y is a fluidity parameter which may depend on time, invariants of the 
rate of deformation, etc., and f denotes a reference value of f that makes the 
expression non-dimensional. To ensure no viscoplastic flow below the yield 
limit we write 


< *(f/f 0 )> 


0 

<M £ A 0 )> 


F < 0 


A 

F > 0 


(3.49) 


A sufficiently general expression for (f> (f/f Q ) is given by a power law 

</> (f/f Q ) = (f/f Q ) n ■ (3.50) 

Various yield criteria and plastic potentials can be introduced depending 
on the material under study. For isotropic materials, for example, there are 
two well-known yield criteria: the Tresca yield criterion assumes that the 

yielding occurs when the greatest difference between any pair of the principal 
stresses cr^, and cr^, reaches a specific value, Y. Usually, Y is taken to be 
the yield stress in uniaxial tension. The von Mises yield criterion for plastic 
yield is 

(<T 1 - a 2 )2 + {a 2 ' ff 3> 2 + ( a 3 - = 2 y2 


or expressing in terms of cm, we have 

^11 ~ °22* + ( 0 22 " U 33' > ^33 ” ^1 1 ^ + 6 ^ a l2 + ^13 + ~ 2 Y (3.51) 


In what follows we shall take a more general yield function that is. defined in 
terms of stress Invariants as 


f^ij) = f(J r J 2 ,J 3 ) 


(3.52) 
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where J^, J" and are the stress invariants, 

J-, = S,-.-, J-= =- ~ S..S.., J- = 4 1 s. • S'- s 

1 ii* 2 2 ij jj.* 3 3 1] ]k ki 

Clearly, von Mises yield criterion, Eq. (3.51), is a special case of (3.47), 

A 

where F is given by 


F <V = y 3 h - YB.F(S..) 


(3.53) 


where Y is the uniaxial yield stress. 


For a special case in which ^ = F and Y is independent of cr.^., we have 
from (3.48) 


d If = y(f/f o ) n 3F/3o-.. = y(f/f Q ) n 3f/3S. j 


and for n = 1, and f given by Eq. (3.53), this reduces to 


d 7? = | (y/f )S.. 

ij 2 w/ o ij 


(3.54) 


After normalizing the yield surface by the von Mises yield stress, i.e., letting 

A 

F = F/Y, a more general form analogous to (3.54) follows from (3.53), namely 


d i? = 


(3.55) 


where 


F = F/Y = 


Ts 


S 

mn mn 

2 2 

j Y 


- 1 


and the plastic potential 


0(F) = 


■ 0 for F < 0 

0(S. .) for F > 0 
ij - 
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Thus, the stress flux for this case becomes upon substituting (3.55) into 
(3.43), while noting = 0, 

S y = 2/i(d’. - y0(F) S../yfT 2 ) (3.56) 

Here, the material parameter y and the plastic potential 0 have to be found 
from experiments. For example, y and 0(F) of 2024-T3 aluminum are readily 
obtained as 

0 = - 1 for 0 < 0 q 

0 +1 

0 = 0 + (F - F ) for 0 > 0 

o 9 o o 


and y = 0.15 yj3 sec \ where 0^ = 10 ^, F q = 9 log(0 Q + 1 ), 
material function which depends on the generalized plastic 


and 9 is a new 


strain e as 
P 


e(e p ) 


0.0003 + 0.0214 e 
P 

- 0.0243 e p 

for 

*^p ^ ® *4 

0.005 


for 

e >0.4 
P 


Here, the generalized plastic strain is 

t 

e ~ f 4 T dP dt 

p J 1 3 mn nm 

0 


3.4.2 Three-Dimensional Failure Criterion for Fiber Composites 

It has been found from experimentation that unidirectional fibrous 
composites exhibit five primary failure stresses, namely, the longitudinal 
tensile,, longitudinal compressive, transversal tensile, transversal com- 
pressive, and intralaminar (in-plane) shearing stresses. In addition, ex- 
perimental observations also show that unidirectional fiber composites under 
combined loading fail at stress levels considerably different from those under 
simple loadings f3o]. Fox generally orthotropic materials, there are four 
more basic material strengths relative to the normal axis, i.e., the normal 
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tensile and compressive, and two normal plane shearing stresses. These 
nine stresses are illustrated in Fig. 3-1. Accordingly, nine simple strengths 
and a combined -stress strength criterion- -are needed to fully describe the 
strength (limit) behavior of the materials. The strength depends not only on 
the properties of the constituents, but also the particular fabrication process 
(such as void content, size and distribution, filament spacing, residual stress, 
etc.). These effects must be included into the formulation of the criterion. 

Various mathematical models have been proposed to describe the strength 
(limit) behavior (cf.. Chapter 2 of Ref. 30 for detailed discussions). In this study, 
we adopt the formulation constructed by Chamis based on a two-level, linear, 
semi- empirical theory [26]. 

The first level is a theory to predict the simple strengths from con- 
stituent properties and fabrication process considerations. Formally, the 
simple strengths can be related to the constituent properties and the fabrica- 
tion process as follows: 


S i = f[(k,d,N,A) f ,k m , (E, y,G,S, e ) ,S B , S ] 

3 * f, m 


(3.57) 


where in (3.57) denotes a certain uniaxial simple strength; (k, d,N, A)^ 

denotes volume content, size, number, and distribution of filaments and voids; 

k denotes the volume content of the matrix; (E, v, G, S, e ') represents the 
111 P f , m 

elastic and strength properties of the filaments and matrix; and and 

denote interface bond strength and residual stress, respectively. The void 

content, the bond strength, and the residual stresses are dependent on the 

filament surface treatment. They also depend on various matrix additives, 


hardeners, temperature, and pressure during fabrication and on the fabri- 
cation method of making the composite. The simple strengths, S^, for various 
cases can be obtained either from the appropriate tests directly, or by evalua- 
tion of the function in (3.57) utilizing the micromechanics structure theory [26]. 
These strengths, together with other material properties, for some typical 
materials are listed in Table 3-1. 
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Table 3-1 

NOMINAL PROPERTIES FOR SOME TYPICAL UNIDIRECTIONAL COMPOSITES 


Properties 

Symbols 

Units 

Graphite 

HM/ERLA4617 

Graphite 
AS/PR5 208 

Kevlar 49/ 
ER LA4617 

S -Glass/ 
ERLA461 

Boron/ 
Avco 5505 

Density 

P 

lb/in^ 

0.056 

0.057 

0.055 

0.077 

0.073 

Elasticity Modulus 

E m 

10 6 psi 

27.5 

18.5 

10.0 ’ 

7.2 

29.2 


E 122’ E i33 

10 6 psi 

i 

1.0 

2.0 

0.6 

2.8 

3.15 

Shear Modulus 

G il2’ G ^13 

10 6 psi 

0.9 

0.6 

0.4 

1.0 

0.78 


* 

G ^23 

10^ psi | 

0.4 

0.4 

0.2 

0.5 

0.60 

Pois son's Ratio 

V JHZ > v m 

— 

0.1 

0.2 

0.4 

0.3 

0.17 


❖ 

^23 

— 

0.4 

0.5 

0.4 

0.5 

0.53 

Thermal Coefficient 

a m 

10" 6 in/in/F 

-0.5 

3.3 

-1.4 

2,1 

3.40 

of Expansion 








a !2Z’ a m 

10" 6 in/in/F 

25.6 

16.2 

. 32.0 

19.9 

18.0 

Uniaxial Failure 
Stresses 

s niT 

lcsi 

122.0 

181.0 

206.0 

257.0 

199.0 . 


s iiic 

ksi 

128.0 

165.0 

28.0 

146.0 

232.0 - 


S i22T' S ^33T 

ksi 

6.1 

8.0 

0.8 

11.0 

8.1 


S i22C’ S jf33C 

ksi 

28.5 

25.0 

9.2 

26.5 

17.9 


S £12 ,S il3 

ksi 

8.9 

13.0 

3.5 

10.2 

9.1 


* s s 

S C3 

ksi 

5.6 

9.1 

2.7 

11.5 

8.9 

Ply Thickness 

S 

*J 

in. 

0.0050 

0.0049 

0.0046 

0.0046 

0.0051 


4 

These values were estimated using composite micromechanics. 



The second level is a theory to predict the onset of failure in terms of 
simple strengths predicted in the first level theory. The governing equation 
for the failure criterion is derived from the following two postulates: 

© At the onset of failure, the distortion energy under simple as well 
as combined loads remains invariant; and 

9 The tensile and compressive properties of the materials are the 
same up to the onset of failure but may be different thereafter. 

The formal derivations are given in detail in Ref. 27, and the resulting equa- 
tion is presented as follows. 


Chamis 1 Failure Criterion for Fiber Composites: Let cr^ denote the 
stress, the uniaxial strength (S > 0}. Then the functional describing the 
three-dimensional, combined- stress failure criterion in fiber composites 
is given by 


r «V S i’ K i> = 1 



(3.58) 

K °1223 g -l33y _ „ a nia °jg33 r K a nig °122B 

i23|37 S £22[3 S £33y il3a7 S in a S i33y n2a $' S nig S ! 22 $ . 


which states that 


F(crg, S^, K^) > 0 No failure 
F((Tg, S^, K^) = 0 Incipient failure 

F(ff^,S^,K^) < 0 Failure 


(3.59) 


Here, are certain coupling coefficients; the subscript l denotes unidirec- 
tional composite (or ply, layer, lamina) properties; the subscripts a, 3, 7 
denote T for tension or C for compression; the subscript S denotes shear; 
and the numerical subscripts 1, 2, 3 refer to an orthogonal coordinate system 
with 1 taken along the fiber direction,2 transverse to it and 3 through the 
thickness. 
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In. the application of Eq. (3.58), it is assumed that the coupling coeffi- 
cients, K^, the stresses, cr..., and the unidirectional strength, S^ are known. 
Following the derivations given in R.ef.2'6, the coupling coefficients. Kg, can 
be expressed in terms of engineering material constants as 


TC 

i23 
K f 13 


• d + 4 ^ i23 - E f22 * ^ ~ ^123^ E ill 

Pill E f22 (2 + V &12 + v n3 i (2 + V £ll * V £Z3^ ^ 
(! -fr 4^ 31 - v n? ) E £33 4- (1 - y jg31 ) E^ 22 

Pi22 E i33 (2 + ^23 + ^21* {2 + V J>32 + V J>31^ ^ 

(1 .* 4v J131 " V $ 32 ) E ill + (1 " y ll2* E i33 
Pi33 E ill (2 + V l3l + V S>32 ) (2 + V £13 + v 112^ ^ 


(3.60) 


where E denotes the elastic modulus, v the Poissons’ ratio. The values of 
E, v, together with the unidirectional strengths are known either from 
measurements or from composite micromechanics; whereas, the stresses 
are obtained by three-dimensional analysis. 

3.5 EQUATION OF STATE 

Since the governing equations contain the pressure (implicitly, in o\ .), 
a constitutive relation must be used relating the pressure to the density and 
the specific internal energy. This particular constitutive relation is well 
known as the equation of state. For high velocity impact of solid bodies, 
two equations of state are well known. One of these is developed by Tillotscn 
[31 ] and the other is developed by Osborne and his associates at Los Alamos 
Scientific Laboratory. 

Till ot son's Equation of State; Tillotson's equaticsi of state (for either 
compression or expansion regions) is given by 

P = ff(e,p) + A[l + B/j, 2 (3.6la) 

for p > 0 with 0 < € < e , and 
^ ‘ o s 
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(3.6lb) 


P = a ep + 


-j + Aft 


-a(l /ri-lY 


(1 + e/e Q rf) 


for p < p Q with. e > where 


p = r) -1) r? = p/p Q , P o = density 


7T ( e, p ) 


ep 


a + 


e/ e Q rj 2 + 1 J 


and a, b, A,B and are parameters (constant) which depend on the material, 
and £ g is the sublimation energy. Values of these parameters for some 
materials are given in Table 3-2. 


Table 3-2 

VALUES OF THE PARAMETERS IN TILLOTSON’S EQUATION 


Parameter 

Aluminum 

Iron 

Copper 

Lead 

p Q (gm/cm' 5 ) 

2.702 

7.86 . 

8.9 

11.34 

a 

0.5 

0.5 

0.5 

0.4 

b 

1.63 

1.5 

1.5 

2.4 

A (mb) 

0.752 

1.28 

1.39 

0.466 

B (mb) 

0.65 

1.05 

1.1 

0.0026 

a 

5.0 

' 5.0 

' 5.0 

10.0 

P 

5.0 

5.0 

5.0 

2.0 

e (mb - cm 3 /gm) 

0.05 

0.095 

0.325 

0.02 • 

° 3 

e s (mb-cm /gm) 

0.03 

0.0244' 

0.0138 

0.0026 


Note that (3.61) is not defined for certain states, for example when p < p Q 
with € < ,e gJ and p > p Q with € > e . Some investigators have used some 
kind of average of the pressures given by (3.6la) and (3.6lb) for these states. 
At this writing Tillotson , s equation is not used in the calculations. 
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Los Alamos Equation of State ; Another equation of state that is widely 
employed in high velocity impact calculations was derived at Los Alamos 
Scientific Laboratory. The equation is given by 


[Afi + 6 E (B + 5EC)]/{6E + <p ), M > 0 

[fiA 1 + 8 E (B q + g Bj + 6EC)]/(SE + (p Q ), M < 0 


(3.62) 


where 

8 = P Q /p, M = P/P Q ~ 1 

A = B = B Q +p{B 1 +pB 2 ) 

C = C^+fiC^, E = pe 

and A, , A^. B , B , . B^, C , C, and (0 are constants that depend on the material. 
1’ 2 o’ 1 2’ o 1 T o 

Values of these parameters for plexiglas, graphite, aluminum, iron, copper 
and lead are given in Table 3-3. 


Table 3-3 

' VALUES OF THE PARAMETERS IN THE LOS ALAMOS EQUATION 
OF STATE FOR VARIOUS MATERIALS 


Parameter 

Plexiglas 

Graphite 

Aluminum 

Iron 

Copper 

Lead 

P o 

1.18 

2.25 

2.702 

7.86 

8.9 

11.34 

A 1 

0.006199 

0.1608 

1.1867 

7.78 

4.9578 

1.4844 


0.015491 

0.1619 

0.763 

31.18 

3.6884 

1.6765 

B 

o 

0.14756 

0.8866 

3.4448 

9.591 

7.4.727 

8.7317 

B i 

0.05619 

0.5140 

1.5451 

15.676 

11.519 

0.96473 

B 2 

0.50504 

1.4377 

0.9643 

4.634 

5.5251 

2.6695 

C 

o 

0.5575 

0.5398 

0.43382 

0.3984 

0.39493 

0.27732 

c i 

0.6151 

0.5960 

0.54873 

0.5306 

0.52883 

0.43079 

'f’o 

0.100 

0.500 

1.5 

9.00 

3.60 

3.300 


The parameters are fitted for gram-centimeter -microsecond system of units. 
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A comparison of the pressures computed from Tillotson’s equation of 
state and Los Alamos equation of state is made for aluminum (p = 2.702 
gm/cm ) at e = 0.5 mb-cm /gm (see Fig. 3-2). The equation of state which 
closely agrees with the experimental data for a given material will be used. 

3.6 INITIAL AND BOUNDARY CONDITIONS 

To complete the description of the high velocity impact problem, we 
must include the initial and boundary conditions, which, in general, are time 
dependent. 

Initial Conditions: At time t = 0 values of all the dependent variables 
(p, v, e , P, S) must be specified at the nodal points of the mesh. It is not 
essential to specify all these quantities at the same set of nodal points. For 
instance, in the analysis of target alone, the velocities, density and specific 
internal energy are specified at one set of points, and pressure is specified 
at a different set of points. 

Boundary Conditions: Depending on the type of the boundary (e.g., rigid 
boundary, free surface, interface, plane of symmetry, etc.), there are different 
kinds of boundary conditions in a problem. At a rigid boundary the normal 
component of the particle velocity must coincide with the normal component 
of the velocity of the rigid boundary. For a fixed (in time) boundary, the 
normal component of the particle velocity must be zero at that boundary. A 
plane of symmetry can be interpreted as a fixed boundary. On the fr.ee sur- 
face, the total stress must vanish. Whereas, the movement of the free 
surface will be included as an integrated part of the coupled Eulerian- 
Lagrangian finite element solution of the governing equations. At an inter- 
face (as well as a contact discontinuity) the total stress and the normal 
component of particle velocity must be continuous, and the density, internal 
energy and the tangential component of particle velocity may be discontinuous 
(jumps). 
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Pressure, P (megabars) 



(P o = 2.702, € =0.5) 
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4. BASIC EQUATIONS IN LAGRANGIAN ZONE 


In this section, the basic equations to be adopted to describe the dynamic 
behavior in the Lagrangian zone are presented. It is clear that, away from the 
impact zone, the disturbance in the target is comparatively small. Thus, the 
conventional structural analysis suffices to deal with the dynamic response in 
this region. 

4.1 THE CHOICES OF COUPLING VARIABLES AND BASIC EQUATIONS 
IN LAGRANGIAN ZONE 

There are various forms of equations governing different types of prob- 
lems in structural- dynamics; yet they consist of different sets of primary 
variables. Hence, the choice of an appropriate set of equations with certain 
primary variables for the Lagrangian zone becomes important in the develop- 
ment of the coupling code. The governing equations and thus their finite ele- 
ment analog to be adopted here must not only be appropriate for tackling the 
dynamic response in the Lagrangian zone for a wide range of problems, but 
the mathematical formulation must also follow the following prerequisitions; 

© In the sense of "coupling," the coupling variables introduced in 
Section 2 must be a part of the primary variables in the governing 
equations, say, e.g., the displacement field, or stress field, hybrid, 
etc. 

© The finite element analog of the equations must be compatible with 
the one in the impact zone, and can be coupled smoothly in the 
mathematical sense as well as the computational sense. 

o The finite element code for its part can either be incorporated into 
the CELFE code as an integrated part of the complete code, or found 
in the existing programs for handling large structural systems such 
as NASTRAN. 

The choice of the coupling variables is not unique in general, but can 
be optimized to reduce the computational efforts. Recall that the primary 
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variables in. the impact zone consist of density, p, velocities v. (or momentums 
V. = p y.), stresses and mesh coordinates x. in a general moving system 
governed by Eq. (3.3). (In the- Eulerian zone contained by the impact zone, the 
last set of variables, x., is of course fixed.) In structural dynamics, on the 
other hand, there is a well developed finite element displacement method, 
finite element force method, mixed models, hybrid models, and free vibra- 
tion analysis, etc., that can be suitably applied for analyzing the response in 
the Lagrangian zone. 


If the force method is used,"e.g., one must choose the stresses (or devia- 
toric stresses) as the coupling variables for the impact zone. In this case, 
the approach complicates the coupling process. The complexities are caused 
not only by the indirect relation on forces and stresses between the Lagrangian 
and the impact zone, but it is also seen in the impact zone that the stresses are 
related to a secondary variable, P, the thermodynamic pressure. This rela- 
tion may lead to a certain difficulty in formulating the coupling procedure 
whenever the stresses .are a part of the coupling variables. 


Thus, in view of the above mentioned prerequisites, the finite element 
displacement method appears to be a better approach. First, there is an 
existing module in NAS TRAN based on this formulation. Second, all the 
primary variables of this formulation, consisting only of displacements and 
velocities, are taken to be the coupling variables for the Lagrangian part in 
the coupling process. In particular, they can be related directly to their 
counterpart in the impact zone. This makes the coupling process rather 
straightforward. Hence, in this study, we will adopt the finite element dis- 
placement approach. 


4.2 ELASTODYNAMIC EQUATIONS FOR LINEAR ANISOTROPIC MATERIALS 


as 


The equations of motion for linear elastic materials are readily obtained 



(4.1) 
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where u. is the displacement along the x- -axis, and F; is the corresponding- 
J 11 

forcing term. 


For fiber composites where the material stiffness coefficients, Ch.^, 
are given in Eqs. (3.40) and (3.41), the equations of motion (Eq. (4.1)) can be 
deduced to the following form: 


a 2 
a u. 


at 


(c 

dx. ^ jjkk 8x k / + 2 


Z 8x, 


/3u. 


[ C jkjk( i 0x k + 9x 


^)] 


+ F . 
1 


(4.2) 


For isotropic materials, the’ equations are reduced to a well-known system 
of wave equations: 

2 2 
a u. ? ' a u. 

p j- = p V u+a+M) „ k -- 4F (4.3) 

ar 3 ■ dx j dx k 3 

where is the Laplacian operator. 


4.3 A GENERAL TRANSIENT ANALYSIS • • 

It should be noted that, in three-dimensional finite element approxi- 
mations of the above equations, three-dimensional solid elements must be 
used. Thus, it turns out to be less efficient for tackling large structures in 
most of the cases under consideration. In order to include the uses of special 
elements, such as bending, membrane, etc., the following system for general 
transient analysis is introduced: 

p 6 + F(S, S;x, t) = G(x, t) (4.4) 

where p(x) may be the density or mass, x = Here, 8 is the de- 

flection vector consisting of lateral as well as rotational components, F is 
a specified vector valued function characterizing the elastic as well as the 
damping effects, and G is a forcing function. 
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When 6 = i.e., 6 does not contain the rotational components. 


1 a I /9u. 3u \ 

fr _ 1 r t L . 

2 ax k . -jki-m iax m 9x-J 


G = F. 

3 

for j, k, £, m. = 1, 2, 3, Eq. (4.1) follows from Eq. (4,4). For general linear 
materials, Eq. (4.4) yields 


p6 + B 6 + KS = G 


(4.5) 


where B is the damping factor, K the spring constant, and 


'v j = 1, 2, 3) 

6 = 3 I 

Qy j = 1,2, 3 J 


Here, 0., j = 1,2,3, are angular displacements. 

J 
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5. A- FINITE ELEMENT ALGORITHM BASED ON THE THEOREM 

OF WEAK SOLUTIONS 


When shocks, compression waves, or other types of discontinuities occur 
in the material, the gradients of the dependent variables in the governing 
equations become very large in the neighborhoods of discontinuities. These 
large gradients lead the effective diffusion terms in the equations to be nega- 
tive in some regions, and cause numerical instabilities in the interaction pro- 
cess during the numerical computations if conventional methods such as 
Galerkin' s, least squares, etc., are used. To overcome these difficulties, 
an alternative finite element formulation is proposed in this section. The 
procedure is utilized primarily for analyzing the dynamic behavior in the 
impact zone. 

The present formulation consists of two primary portions. First, the 
finite element formulation is constructed based on the theorem of weak solu- 
tions, so that the jump conditions can be satisfied automatically to take care 
of shock propagations. Second, a generalized two-step, time -splitting, shock 
capturing scheme is developed and implemented. The scheme so construed 
has a capability to remedy the spurious oscillations arising due to numerical 
instabilities. 

In Section 5.1, a general discussion on the theorem of weak solutions 
will be presented without proof. The governing equations in conservation 
form and the finite element analog of these equations, are presented in the 
subsequent subsection. The formulation of a two-step scheme is discussed 
in Section 5.3. 
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5.1 ON THE THEOREM OF WEAK SOLUTIONS 


The results from the theorem of weak solutions have been frequently 
applied in developing finite difference schemes for solving first order, non- 
linear hyperbolic equations. In most of these schemes, the advantages of the 
resulting forms deduced from the theorem are not fully utilized, since they are 
always in the integral form. This fact, however, can be implemented in a 
simple way into a finite element scheme, and the advantages can be fully 
utilized. 


Various works on finite difference have indicated that, when the results 
of the theorem of weak solutions are to be used in the numerical scheme, the 
governing equations need to be converted into a conservation form {cf. Ritchmyer 
and Morton [32] ). In what follows, therefore, we restrict ourselves in consider- 
ing the system of hyperbolic equations in conservation form. It is obvious that 
all the conservation equations can be written in this form. The constitutive 
equations, however, are not in the conservation form hut can he recast into 
such form as shown in subsection 5,2. 


5.1.1 Theorem of Weak Solutions for First Order Quasilinear Equations 


Consider a Cauchy problem of a system of quasilinear hyperbolic equa- 
tions 


M , _ 

9t + 9^ ~ 


k = 1,2,3 


(5.1) 


on a cylinder R = (x) x T (t), where x = (x-^x^Xg), F k = F k {x,t, <£), and G = 

G (x, t, <£). A class of piecewise smooth and piecewise continuous, vector val- 
ued functions $ in R for t > 0 are called the weak solutions of (5,1) if the fol- 
lowing relation is satisfied: 



df2dt 


-If 


i Gdfidt 


R 


(5.2) 
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where t, e C with compact support on R, i.e., is any continuously 
differentiable function that vanishes on 8R, the boundary of JR. 

Lax [33], Oleinik [34] et al., have shown that, if the following conditions 
are satisfied: 


1. For each k = 1,2,3,F has a continuous partial derivative 

.k 


with respect to <j) , and- 


8F 

~df 


is bounded for (x, t) CR; 


q 2 ^k ^2 F k 

2. For bounded <j), are continuous, 


8^ F^ 
and > 0; 

80 - 


8 * 8 V 8 # 2 


3. The function G (x, t, 0 ) has a continuous partial derivative 
with respect to 0, 


the generalized solution of Eq. (5,1) with a piecewise continuous initial condition 
is unique, and satisfies Eq. (5.2). 

The proof of the theorem can be found in the aforementioned references, 
and will not be presented here. However, we will discuss in what follows 
some immediate consequences of the theorem directly related to our proposed 
formulation. 

5.1.2 Jump Condition 

Following the Green's theorem, Eq. (5.2) can be deduced to a form 
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~ fjd t * X ’ 


t) $ (x, t)]' df2 dt + 


R 


fl 

8Q xT 


n k dS dt 


-// 


dft dt 


(5.3) 


R 


where n k is the k component of the outward unit normal vectqr on 9ft, the 
boundary of ft. 

Observe now the integrals 

I v = ff [Z(x,t) tf(x,t)] dftdt (5.4) 

R 

for some t^eT, xeft; and 


=// 

9ft xT 


% ‘(x, t) F 


(x, t) n, dS dt 


(5.5) 


for all xeSft, and teT. Since £, (x, t) is piecewise continuous function with 
compact support in R, it is obvious that I — 0 if there is no discontinuity in 

' • S lr 

R. Otherwise, we may deduce without difficulty the jump of the function F 
across the discontinuity, i.e., 

. i[F k ln k dSdt (5.6) 

8Q xT 
s 

denotes the jump of 


where d ft g is the surface of the shock layer, pr^Jj 
the component of F, across the shock. 
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For the integral I , if there is no discontinuity on R, Eq. (5.4) can be 
reduced to 


I 

v 


■/ 


£ (x, t) (j) (x,.t) 


-/ 


t, (x, 0) (x, 0)dft 


Q 


a 


Across a shock, however, the Leibnitz's rule yields 


/ 


I v = / ?. (x, t) (j) (x,t) d Cl - g £(x, 0) <j) (x, 0) dQ 


■/ 




n 


+ 

« 

T 


■y ^ ^ (xj,d y c s (*n,e) ds(ii) d§ 


9Q (t) 




cue) ds(ri) de 

b 


' T 


on (t) 

s ' 


where the subscripts 1, 2 denote the values upstream and downstream of 
the shock, respectively. Here, the speed of sound is defined as 

C = 
s 


9S 

at n k 


The last two terms in the above integral are no more than the jump of < f>, i.e. 


I v ~ J £ (x, t) $ (x,t) d Q ~ J' t, (x, 0) <j) (x, 0) d& 




Q 
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+ 


(5.7) 



dn xt 
s 


C (x.t) C s (x,t) ['<£] dsdt 


Since the entropy condition requires that 


C s m - |[F k ! ^ = 0 (5.8) 

The substitution of (5.6) and (5.7) into (5.3) then yields 


ff<* lr+ I ' k ^) dndt 

R 



t (x, 0) <b (x, 0) d£2 - 


If 


R 


l, G dfidt 


(5.9) 


This consequence implies that (5.9) satisfies the jump condition (5.8) auto- 
matically. 


5/2 FINITE ELEMENT ANALOGUE OF WEAK SOLUTIONS 
Recall that we are considering a system of equations 




(5.1) 


1c j 

where F = F (x,t, <P), k = 1, 2, 3 and G = G (x, t, <£). In this section, we shall 
formulate the finite element analog of equations in the form of (5.1). As afore- 
mentioned, the present formulation is developed primarily to deal with the be- 
havior in the impact zone. Thus, we have to write the conservation equations 
and the constitutive equations for the impact zone in the above form. 
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5.2.1 Governing Equations of Conservation Form in Impact Zone 

It is obvious that the conservation equations (3.11) through (3,13) are in 
the conservation form. The constitutive equations (3.39), however, are not in 
the form of (5.1). Therefore, we recast the constitutive equations as 


as.. a<a S..) 
+ k « 


at 


dx, 


— = 0S. . + C 


- S. co . + S .co. 


xj ljkf kf lm mj mj lm 


(5.10) 


where © = div fi. Here, an additional term©S.. appears on the right-hand side, 
which is considered as the forcing term. 


In view of (3.11) through (3.13) and (5.10), we may replace the forcing 
function G(x, t, 0) in (5.1) by A <j) + G(x, t) without causing any confusion. Thus, 
instead of (5.1), we have 

k 

If- +H- = +G(x,t) (5.11.) 

k 

For the impact zone, the vector valued functions in (5.11) represent, respec- 
tively, the following physical quantities: 


= < 


v,; j = 1 , 2 , 3 

J 

E 

^ ; i, j — 1,2,3 


F = 


fik P 

°k V j " CT jk ; J = 2 » 3 

a k E “ a ik v i 
f2k ^ >*->j = l>^»3 


; k,i = 1,2,3 
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3x, 


“ v t) 


k 

a 


A<£ = 


1 v j a% Pk - V' i “ ! - 2 - 3 


E air »k- V 

k 

3Q. 

' s ij aS? ; ^ = 1 > 2 > 3 


; k = 1,2,3 


G = byj = 1 ’ 2 ’ 3 

F i v i 

C ijki d ki " S im “mj + S mj W im J if 3 = *• 2 ’ 3 


l;k,i = 1,2,3 (5.12) 


5.2.2 Finite Element Formulation 

Suppose that the discrete approximation is constructed by appropriate 
interpolation functions in each element, and take £(x, t) defined in (5.9) to be 
the weighting function with compact support in the cylinder AR = Q e x T. Here, 
is the volume of an element e, T = [o, AtJ and At is the time step in the time- 
split scheme. In order to utilize integral relation (5.9) to minimize the error, 
we have to approximate the (weighting) function £ by the shape function at the 
previous time step, i.e., 

i T = (1 -^j)N r , r = l,2,...m (5.13) 

where m is the number of nodes in the element e. Also, approximate the 
solution to (5.11) as 

* = N s [ (1 ‘ + Zt^i n+1) } s = 1 »•*•», m (5.14) 
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With these expressions, we obtain a matrix equation in the form of 


(KJ - f ■ ((*„] * “ Ml* 


(n) 

s 


+ t{Cj (5.15) 


where 


rs 


/ 

a. 


N N dO 
r s 


rs 


■/[ 

■f\ 


e 

9N 

+AN r 


, 3N : 

H &r + GN r 


N dO 
s 


da 


(5.16) 


Here, the vector H = F - -6Q . 


5.2.3 Finite Element Analogs of Equation (5.11) in Eulerian and Lagrangian 
Modes 

Eulerian Representation; When = v ] s -» k = 1> 2, 3, the conservation 
form (5.11) becomes 


M + = G 

at ax, 

k 


(5.17) 


with 


= < 


V,; j' = 1,2,3 

J 

E 

. S. i, j ~ 1, 2, 3 
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F = { 


V k P 

V k V j - a jk ; j= 1,2,3 
V k E '- ^k V i 
v k S ij ; J = 1> 2 > 3 


; k, t = 1 , 2-, 3 


G = 


F.; j = 1, 2, 3 

F v 
JL 2 

01 

. C . . , . d. - S . 60 . + S ,to. ; i, i = 1, 2, 3 

l ljki ki im raj mj im J . 


; k, i = 1, 2, 3 


The corresponding finite element analog of (5.17), (5.18) still retains the 
similar form as (5.15), but with different coefficient matrices, i.e.. 


(M -t? [»„])(*ri ■ ([*„]*¥ mi 


{^i n) j + At jc 


where 


A 


rs 


/ N r N s d Sl 




B 


rs 


■/ 


9N 


k 9x k s 


N dfl 


- /[* k < 
n„ 1 k 


+ N r G 


dO 


(! 


Lagrangian Representation: When£l k — 0, k = 1, 2, 3, the functions F 
and A are reduced to the following forms: 


F = 


-<7 jk’ 3 ” 1, 2, 3 


-Ob, v„ 
ik i 


V ; k, i = 1, 2, 3 


(5.18) 


(5.15) 


>.19) 

,k 
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A<£ 


-P 


3x, 


Ov, 


■V. 


j ax. 


■; j = 1,2, 3 


' -E 


3x, 


0; i, j = 1, 2, 3 


; k = 


1,2,3 


(5.20) 


For this case, the corresponding finite element analog does not retain the 
original form. Instead, it yields the classical form of explicit scheme: 


with 


Me +1) j 


• rs 


/ N N d& 
r s 




c - ■/ 1 - 


k 9N r 

F -s-i + N G 
3x k r 


dQ 


(5.21) 


(5.22) 


5.2.4 Remarks 

Instead of using Eq. (5.13), we may approximate the weighting function by 
£, r = t/At N , r = 1,2,. . '. , m, as well. In this case, however, the integral formula 
(5.3) instead of (5.9) has to be used. Mathematically, these two formulations could 
achieve the same accuracy and possess a similar stability character. Neverthe- 
less, the surface integrals appearing in (5.3) would affect the numerical computa- 
tions due to the propagative nature of the hyperbolic system. 

It is important to note that neither the use of - t/At N r nor £ = (1 - t/At) N r 
can always represent accurately the class of test functions having compact support 
on R if the conventional shape functions are used. These known classes of shape 
functions may not vanish or be equal to the same constant on all boundaries ' 
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of the cylinder, R, although we have presumed that they do as required by the 
theorem of weak solutions. The error arising from this aspect, nevertheless, 
can be reduced easily by treating the boundary conditions with care during the 
numerical computations. Meanwhile, constructing a new class of shape func- 
tions does not seem to be practical due to the difficulties involved, and the 
serendipity elements will be used instead. 

5.3 ISOPARAMETRIC ELEMENTS AND NUMERICAL INTEGRATION 

As is seen, the use of finite elements discretizes a continuum problem 
and establishes a system of algorithm equations, whose coefficients are ex- 
pressed in terms of the products of element shape functions. The choice of 
element type in finite element analysis is usually dictated by considerations 
of accuracy, computational efficiency, and the specific problem under study. 
As remarked in the previous subsection, the serendipity family will be used 
in the present study. 

For three-dimensional finite element analysis, the "serendipity" ele- 
ments with isoparametric formulation are superior to other solid elements 
with respect to the above considerations especially in the present problem 
where an adequate representation of the geometry is essential. The term 
"isoparametric" means that the same shape functions are used to' define both 
the geometry and the unknown function. In the present study, use will be 
made of the linear, quadratic, and cubic elements of the "serendipity" family. 
After performing certain transformation (mapping), these brick-type ele- 
ments will deform to yield curved surfaces. The following is a description 
of these elements. 

5.3.1 Serendipity Elements 

In the "serendipity" brick element, most of the nodes are located on 
external edges. In fact, the linear, quadratic, and cubic elements contain no 
internal nodes at all. The corresponding shape functions are listed on the 
following page: 
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Let 

*o = ™i and S 0 =Ki 

where (£ } rj , £) are the local coordinates and (£^, rj ^ j^) denote the coordinates 
of nodal points of 'the cube (see Fig. 5-1). 

Linear Element (8 nodes) 

The shape functions in this case are defined by 

N. = 1 (1 + £ ) (1 + IJ 0 ) (1 + i Q ), i = 1, 2, ,8 

8 

Quadratic Element (20 nodes) 

For corner nodes 

N i = f< 1 + M 1 + V( 1 + h><£o + n 0 +t 0 -2) 

Typical mid -side node 

t i = 0, rj i = + 1, C i = + 1 

Cubic Element (32 nodes) 

For corner -nodes 

N i = M {1+ § o )(1 + V (1 + ^ 0 } [ 9{|2+ 1l2+ ^ 2) " 19 
Typical side node 

f ^ = ± ^ + 1 

N i = m (1 -e 2 )(i + 9e 0 )(i + T lo )(i + ^ o ) 
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5.3.2 Numerical Integration. 

In finite element analysis, the matrices defining element properties, 
e.g., stiffness, etc., must be found. These will be of the form 


I 


Iff 


[G(x, y, z)] dxdydz 


(5.23) 


in which the expression G depends on the equation being solved and the shape 
functions N^ and/or their derivatives with respect to the global coordinate 
system (x, y, z). As is seen earlier, the shape functions are written in the 
local coordinates (£, q, £,), therefore certain transformations must be per- 
formed so that Eq;(5.23) can be evaluated in the local coordinate system. 

First, the expression G(x, y, z), which involves shape functions and 
their derivatives in the global coordinate system (x, y, z) must be transformed 
into those in the local coordinate system {§, rj,^}. When isoparametric formu- 
lation is used, the relationship between global coordinates and local coordinates 
is defined by 

X = N. X. 

y = N. (£ , Ti , & ) y. (sum on i) (5.24) 

and 

z = N. (£,ti,U z. 

in which are the shape functions as defined before and x^, y^, z^ represent 
the global nodal coordinates. By this transformation, the originally right 
prism in the (£,r|,£) space will become distorted in the (x, y, z) space. Now, 
since the shape functions are defined locally in finite element analysis, no 
transformation is necessary. However, the derivatives of shape functions 
with respect to (x, y, z) must be transformed into the local coordinate system. 
This can be done as follows. 
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Using the chain rule, one has 


3 x 3 y 3 z 

3 3 1 , 0 £ 


where |jj is the Jacobian matrix. Therefore, one obtains 


The Jacobian matrix ^jj and its inverse, in turn, can be determined from 
(5.25). More explicitly, the Jacobian matrix becomes 


l J J ■ 


l Z 72 z 2 


and its inverse can be determined subsequently. 
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Secondly, a transformation for the volume element must also be done, 
that is 


dxdydz = det J d§dr]d£ (5.27) 


By combining (5.26) and (5.27) > one finally obtains, in place of (5.24), 
the following integral form ' 



G(g, ri, u 


d£ dp, d£ 


(5.28) 


Thus the integration is carried out within the right prism and not in the 
complicated distorted shape. 


While the limits of the integration are simple’ in (5.28) , unfort un ately the 
explicit form of (G) is not. Therefore numerical integration usually has to be 
resorted to. Essentially, (5.28) is approximated by the following form; 


I = 


n n n 

rzE 

i = l j =1 m = 1 


w. w. w 
i j m 


e(g.. n r C m )] (5.29, 


In the above, w., w. and w are' the weighting coefficients with G evaluated 

i j m b 

at the point (£ ., £,^). Herein, for simplicity, the number of integrating 

points in each direction was assumed to be the same. The numerical scheme 
presently used is the Gaussian quadrature, of which the points for evaluating 
G and hence the corresponding weights, are preselected to yield higher accu- 
racy with a fixed number of integration points. Of course, within the iso- 
parametric family, exact numerical Gaussian quadratures can be obtained 
if a sufficient number of points is used. 


The surface integrals can be expressed in a similar form 


I = JJ [H(x,y,z)JdS 


(5.30) 
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in which H(x, y, z) is obtained from the shape functions and/or their derivatives 
with respect to the global coordinate system (x, y, z) and S is a curved surface 
in space. Since the integrand is us.ually very complicated, expression (5.30) 
will be evaluated by numerical integration. 

First, the expression H(x, y, z) must be expressed in terms of the local 
coordinates (£, p, £), analogous to the volume integrals, and to be evaluated 
on the appropriate surface. Second, a transformation for the area element 
must be done so that the integration is .performed on the surface defined by 
two of the local coordinates. For instance, the area element on a surface 
where £ is constant can be written as 


dS = detff] d£ dr) (5.31) 

in which [J] is the modified Jacobian matrix evaluated on the surface con- 
sidered. The matrix, in turn, is defined as 



(5.32) 


with 



The corresponding surface integral, througn the above transformations, 
thus becomes 


I 


If fH(g.n)] d£d n 

-1 -1 


(5.33) 
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Finally by using Gaussian quadrature, the surface integral is approximated 
by the following summation 

n n ’ . 

1 = £ £ W i w j (H^i’TJj)] (5.34) 

i=l j=l 

in which H is evaluated at the Gaussian points (£., n .), w.,w. are the corre- 
sponding weighting coefficients, and n denotes the total number of Gaussian 
points to be used in each direction. 

5.4 GENERALIZED TWO-STEP, TIME-SPLITTING SCHEME 

It is well known that the finite element method with conventional assembly 
techniques is equivalent to the centered finite difference scheme. When the 
method is applied to solve flow problems containing shock waves, numerical 
instabilities generally arise primarily due to the lack of dissipative terms. 

In the finite difference approach, this difficulty can be overcome by either one, 
or a combination, of the following schemes; 

• Introducing artificial viscosity 

a Replacing the center difference by a noncocentric 
difference scheme, 

• Utilizing the two- step. finite-difference scheme, or 

e Introducing Lax-Wendroff' s second order correction, 
etc. 

All of the above variations except the second have been adopted in the current 
finite element codes to solve the impact problem. It was found that the last two 
approaches yielded promising results. However, due to the excessive compu- 
tation time needed with the Lax-Wendroff second order scheme, we have 
decided to adopt the two-step, time -splitting procedure. 

Consider a non-lineal matrix equation 

<[ A rs ] - A* e I B r s!> M n+1) } * <I A rs ] +At 11 - 9) 

+ At{c r } (5.35)' 
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where 9 is a non-negative factor, [A ], [B ] are square matrices, ]C { is 
a column matrix, and denotes the solution of $ at the n step.' In the 
present cas.e,. the, factor 9 = 1-/3, -the matrices. are 


and 


r s 


= / N N dfi, 

J r s 


B 


/ 


n 


a- 


9N 


% 9x, + A N i 


r, ,s = 1, 3, ...,f 


N dfl 
s 


and the vector 


/* ( 3N 

= / 

d 1 k 


H -f N ? G 1 


where = F^ - v^.0 . 

Following the general approach of two-step procedure (see Richtenyer 
and Morton [32]), the solution of (S.35-) can be obtained by the following steps: 


([Aj^a At e[B< n s ’])j*< n+a) j= ([A rs ] + aAt(l-8)[B<“>| j^’J + aAt jc[ n > j 


j ‘ j' 


([A rs ]-At8 [>«>]] . ([A rs ] + At(l-0) [>«’]) j* 


(5.36) 


* n) |+At lc< n+1) l 

I ( 1 


where 


B (n+1) = cxB <n) + (l-c t )B (n+s > 

■W P * ' VC 


rs 


rs 


rs 


~(n+l) , aC (n) ( 1 . a) c< n +a) 


■1 

Here a is a relaxation factor, and the constant a may be an integer or a frac 
tion multiplier of the time step. As will be seen in Section 7, the numerical 
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computations become more stable as the parameter a increases. The mathe- 
matical aspect of the present formulation is yet to be studied in-depth. In 
particular, the relation between a and the ratio of the time step size and the 
element size must be determined in order to obtain an optimal accuracy for 
the general impact problems. 
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6. COUPLED EULE RIAN - LAGR ANGIAN 
FINITE ELEMENT (CELFE) CODE 


The complete CELFE program to be constructed follows the require- 
ments discussed in Section 2, i.e.^ it contains the capabilities of 

• Simulating the transient state during the impact process and right 
after the process is completed 

• Computing the dynamic response of the entire structure for specified 
times after the impact process is completed, and 

• Detecting various criteria during the impact process, e.g., rebound- 
ing, sliding, penetration, together with tracing the movement of the 
free surface, failure front, etc. 

In the present study, however, we are mainly concerned with the development 
of the transient analysis modules. This aspect will be emphasized in the 
following discussions. 

6.1 BASIC IDEAS AND FOR MUL ATION 

In this section, a preliminary discussion of the procedure for transient 
analysis is given. We shall regard the procedure as the principal CELFE 
procedure or simply the CELFE procedure. 

The CELFE procedure can be divided into three parts: the first part 
contains a finite element module based on the algorithm formulated in Section 
5 to simulate the dynamic behavior in the impact zone; the second part is a 
classical finite element procedure for analyzing the dynamic response in the 
Lagrangian zone; and the third part consists of a coupling (interfacing) procedure 
to integrate the above two separate parts to form a complete system, i.e., the 
complete CELFE system. 
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The complete CELFE system to be developed in this section consists of 
two options, namely, an in-core CELFE code and a coupled CELFE-NASTRAN 
code. The in-core option is built as a substructure of the coupled option to 
allow more flexibility in utilizing the computer code. 

6.1.1 Definitions and Preliminary Discus sions 

Recall that the entire structure has been partitioned into an impact zone 
and a structural (Lagrangian) zone. The impact zone consists of an Eulerian 
zone (E zone) and a transition zone (E-L zone). 

Based on the above discussions, the Lagrangian zone (L zone) can be 
divided into an in-core and a NASTRAN zone according to the options of the 
in-core CELFE code and the coupled CELFE-NASTRAN code. The zonings 
can be denoted as L c zone and L^ zone, respectively. Since the in-core option 
is contained in the coupled option, the Lagrangian zone (L zone) may be com- 
posed solely of L c when the in-core option is used; or composed solely of L^ 
or mixed L c and L^ when the coupled option is used. 

In the entire impact zone, together with the L c zone, three-dimensional 
solid elements have to be used. In the remaining portion of the structure, i.e., 
the L-j^j. zone, any appropriate elements {depending on the nature of the struc- 
ture) can be applied. A typical arrangement of the elements is illustrated in 
Fig. 6-1. Here, once again we stress that; 

• In the impact zone (i.e., E zone and E-L zone), the dynamic be- 
havior will be handled by the algorithm developed in Section 5 
using the hydro-elasto-viscoplastic ‘model represented by a gen- 
eral moving coordinate system. The primary variables in this 
zone consist of the density, p, momentum, pv., j = 1, 2, 3, total 

energy, pe, and deviatoric stresses, S„, i, j = 1, 2, 3; the remaining 

variables, i.e., the pressure P, stresses ct„, and the displacements 

Uj , etc., belong to the secondary variable. 

• In the Lagrangian zone composed of the L c and L^ zones, tradi- 
tional structural dynamic analyses will be applied. In this study, 
the displacement method is adopted, and thus the displacements 
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Interface Projectile 



Fig. 6-1 - Typical Finite Element Sketch in Global Analysis of High Velocity Impact 







(both the lateral and rotational displacements) are chosen to be 
the primary variables. 

• Since the governing equations,, coordinate .systems- and- types -of 
elements in the various zones are different in nature, the coupling 
procedure must provide two mechanisms, namely, an in-core 
procedure to couple the E-L and L c zones (or zone in case of 

empty L^); and an interfacing procedure to couple the main CELFE 

program with NASTRAN. The coupling variables for the present 
procedure contain lateral displacement and velocity fields. 


The associate variables, types of elements, and related aspects concerning 
the above discussions for various zones are summarized in Table 6-1. 


6.1.2 Governing Equations 

For the transient state, i.e., during the penetration process and right 
after the process is completed, the dynamic behavior in the impact zone and 
in the Lagrangian zone are governed by different sets of equations as dis- 
cussed in Sections 3 and 4. The equations are recalled as follows: 


For the impact zone, the equations are represented in a general moving 
coordinate system by the conservation form 


where 


If = M +G(xit) 


V ; j = 1, 2, 3 
E 


^.1.1 -1.2.3 


( 6 . 1 ) 


F = 


f« k P 

Q k V j - CT jk ; j = 1 ’ 2 ’ 3 

n k E -°ik r i 
*2 k S"» i> j = 1> 2, 3 


\ ; k,i = 1,2, 3 
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Classification 
of Zones 

Impact Zone 

Lagrangian Zone 

E-Zone 

E-L Zone 

Lr Zone 
c 

Zone 

Eulerian Zone 
(Empty at t= 0) 

Transition Zone 

In-Core Lagrangian 
Zone (May be Empty) 

' NAS TRAN Zone 
(May be Empty) 

Primary Variables i 

p,pv.,pe,S..; i, j = 1, 2, 3 

J 

u., 0.; j = 1, 2, 3 

J w 

Secondary 

Variables 

u j' p ’ a iy l * j = z ‘ 3 

“r V 3 - 2 * 3 ( "j ' V 

Type of Elements 

3-D Isoparametric Elements 

,3-D Isoparametric 
Elements 

Any Compatible 
Elements 

. 



Coupling Variables 


v., u.; j = 1, 2, 3 

J J 


























ax, ^ k " V k } 


V. 


*k 

a 


(f2. - v, ) 


A<j> = 


j 9x. v “k k 


9x, (il k “ v k J 


^k 

an, 


; k = 1,2, 3 




ij ax 


- ; i, j = 1, 2, 3 


G = { 


F.; j = 1, 2, 3 

J 

F - v 
■i i 

C.. T . df - S. M . + S .00. ; i, j = 1, 2, 3 

ljki ki lm mj mj lm 


L ; k.i = 1,2,3 (6.2) 


For the Lagrangian zone, the small deformation is assumed in this study, 
and the governing equations are given as 


p<5 + BS + KS = G 


(6.3) 


In the L c zone where three-dimensional isoparametric elements are used (cf. 
Table 6-1), we have 


A 

P 


= {u.;j- 1.2. 3} 


pv. = V.; j = 1,2,3 


B = 0 


K = 


1 9 


av ) * i = i> ^ 


2 9x k 'jkjk 9x k ; u ki 


; k,i = 1,2,3 


G = : i~" (C 


au 


k ) ) +F.; j = 1,2,3 ; k = 1,2,3 (6.4) 


ax^ jkjk ax k ; 2 8x k ' jkjk 9Xj ' j 1 
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where denotes the Kronecker delta, and the damping term has been included 
in the forcing term G. In the L^ zone, on the other hand. 


ju.; j = 1, 2, 3 

5 - I ^ 

e • j =.i,2, 3 

" J 


(6.5) 


and p, B, K and G are given according to the nature of the problem. 


6.1,3 Coupling Procedure 


As mentioned previously, the coupling variables to be used in the CELFE 

code are u. and v., j = 1,2, 3. This set of variables provides a straightforward 

1 3 

coupling procedure, i.e., the impact zone may be considered as a substructure 

of the Lagrangian zone where the global finite element matrix equation for u. 

and v. contain a submatrix equation representing the impact zone. This coup- 
3 

ling procedure is called here the direct method. One may also use the sub- 
structural approach developed by Przemieniecki (Ref. 35), which we may call 
an indirect method. 


It will be seen that the direct method is simpler, and the in-core coupling 
procedure proceeds automatically in the CELFE program following the transi- 
tion of the relative velocity SI ^-from v^_ to 0 for k = 1, 2, 3. There is a drawback, 
however, in that its in-core capability is limited by the computer storage. Since 
our objective is to couple CELFE with NASTRAN in the first place, this draw- 
back is not important. 


Although the indirect method may extend the capability of the in-core 
CELFE option (in exchange for increasing the computation time due to the 
complexity in the coupling process), the size of the problems that can be 
handled will still be limited compared to NASTRAN. Thus, this approach 
turns out to be impractical for the present study. 

In view of the above discussions, the direct method is chosen here for 
the coupling procedure of the CELFE code. By definition, the displacement 


6-7 ' 



field is given as: 


( 6 . 6 ) 


u = v., j = 1, 2, 3 

J J 


In the impact- zone-, however, the- velocity field is represented by a general 

moving coordinate system, u. / x. and follows the identity (Eq. (3.3)). From 

3 J 

the expression (Eq. (3.9)), Eq. (6.6) yields a conservation form 





A u. + u., 

J J 


j = 1, 2, 3 


where 


A = 




is the divergence of the relative velocity. 


(6.7) 


Across the impact zone, the relative velocity k = 1, 2, 3, governed 
by Eq. (6.7), will change from = v^ in the E zone in such a manner that, 
in E-L zone, 

n = N 

b 


(1 - 


~)fi (n) 
At ; “s 


+ 4n (n+i) 

At ~s 


( 6 . 8 ) 


for s = 1, . . . , m. At the nodes on the failure front, £2^ begins to vanish and 
remains at zero up to the edge of the impact zone, i.e., the interface between 
the impact zone and L c (or L) zone. Then, the governing equations (6.7) are 
replaced by (6.3). throughout the Lagrangian zone. 

Noting that Eq. (6.8) can be verified by the definition Eq. (3.3), the 
relation ensures that at any instance the mesh in the E-L zone will remain 
simply connected. In addition, since the finite element approximation is 
taken as the average of the solution-field elementwise, Eq. (6.8) ensures 
that the solutions in the E-L zone will maintain the same accuracy as in the 
remaining part without adjusting the manipulations explicitly when intrusion 
of the Lagrangian mesh into the originally Euler ian mesh occur. This sig- 
nificantly simplifies the program for the coupling procedure compared to 
finite difference techniques. 
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6.2 FINITE ELEMENT FORMULATION FOR TRANSIENT ANALYSIS 

The principal function of the CELFE code is to simulate the dynamic 
behavior during the penetration process and right after the process is com- 
pleted. The finite element formulations to be discussed in the following are 
directed mainly for this principal function. 


6.2.1 Finite Element Module for Impact Zone 

The stated module consists of a procedure to solve 'the system (Eq. (6.1)) 
with, listed primary variables (cf. Table 6-1). The finite element analog of 
Eq. (6.1) is readily obtained, as shown in Section 5, based on the theorem of 
weak solutions, i.e., 



I A rs] 


2 At 






+ At 


Ki 


(b.9) 


where 
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- fl, <j>, and the primary variables 
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P 

V-; j = 1,2,3 
3 
E 

S..; i,j = 1,2,3 
. U 


( 6 . 10 ) 
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The system (Eq. (6.9)) is solved independently (in the impact .zone) 
from the remaining zone. Its solution provides a finite element approxima- 
tion of the nonlinear behavior in the vicinity of the impact point due to large 

deformation, in particular, the failure. The momentums, V-, j = 1, 2, 3, in 

3 

addition, will serve as "known 11 functions, dividing by p, for the coupling pro- 
cedure characterizing the velocity field in the impact substructure. 

6.2.2 Finite Element Analogs for Displacements in Various Zones 

Since the coupling procedure is formulated by the direct method intro- 
duced in Section 6.1.4, a global system of finite element equations is required 
for the displacement field covering the entire structure. This finite element 
system is constructed by summing up the subsystems deduced from Eq. (6.7) 
and Eq. (6.3) for the impact zone and Lagrangian zone, respectively. 


Impact Zone: The finite element analog of Eq. (6.7) can be obtained 
following the weak solution formulation discussed in Section 5, namely. 



<C A rJ + t B r J> K 


(n)j 

s 


+ At C 


( 6 . 11 ) 


Here, 



C 

r 



v. N di7 
J r 


( 6 . 12 ) 


and the velocity field v. in is 
the previous subsection. 


determined by the formulation presented in 


6-10 



Lagrangian Zone ; Since the form (Eq. (6.3)) depends on the model used 
for a specific problem, it can be approximated using different types of ele- 
ments with different forms associated with the model. Numerous forms have 
been developed and are available in the literature. Symbolically, the form is 
represented by 

[M] [V] + [B] [ 5 }+ [k] jej = [g| . (6.13) 


The integration with respect to time in Eq. (6.13) can be obtained either by a 
finite difference or finite element scheme. In NASTRAN, and thus in the L^ 
zone of the present analysis, the finite difference procedure is used (cf. Ref. 
36). Whereas in the L c zone, where three-dimensional solid elements are 
utilized, the finite element in time formulated by Zienkiewicz (Ref. 37, Section 
16.5, pp. 335-342) is preferred. 

Let p represent either finite element or finite difference discretization 
of d/dt. Then Eq. (6.13) yields 

j[p 2 M] + [ P B] + [K]j jsj = .{gJ (6.14) 


6.2.3 Finite Element Formulation of Coupling Process 

As mentioned previously, when the direct method is adopted to couple 
the Eulerian and Lagrangian modes, it is necessary to formulate a global 
system of matrix equations for the displacement field in such a way that: 

s The system covers both the impact zone and the Lagrangian zone, and 

• The system is constructed by summing up the finite element analog 

(Eq. (6.11)) when the element belongs to the impact zone and Eq. 

(6.14) when the element belongs to the Lagrangian zone. 

Let a superscript, e, denote the element value of a given physical quan- 
tity, and subscripts E-L and L denote the corresponding quantity defined in 
the impact zone and Lagrangian zone, respectively. Let also 

[Kj] - [Aj - f [BJ 
J E-L 
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! s 4-l = ([A - ] +iL ^ [B - ]) 8:1 + At K! 


for j = 1, 2, 3, where [a ], [B ], and |C j are defined in Eq. (6.12). Then, 

* S TS ( JT j 

Eq. (6.11) can be written as 


[ K j] h| (e) = [Gj! ( * 

i E-L { J, E-L 1 


e) 

’E-L 


(6.15) 


Similarly, we may write Eq. (6.14)’ in the following form 


j[p 2 M E+. [p b ] l + Ml| H 6) = ia!<e) 


(6.16) 


Summing Eqs. (6.15) and (6.16) over e = 1, 2, . . . ,M, yields 


j[p 2 M] + [pB] + [K]j jej = ja 


(6.17) 


where M is the total number of elements in the entire structure, including 
the projectile. 

Noting that the actual coupling variables consist of lateral displacements 

and velocities, u. and v., j = 1, 2, 3, in the impact zone the lateral velocity 
J J 

(momentum) field is categorized as primary variables, and is obtained follow- 
ing the procedure given in Section 6.2.1. Whereas, the lateral displacement 
field belongs to the secondary. In the Lagrangian zone, on the other hand, 

Eqs. (6.3), (6.4) and (6.5) indicate that, in general, there are six degrees of 
freedom (i.e., six primary variables) composed of three lateral and three 
angular displacements. The velocities, in turn, become secondary. 

Therefore, when we sum Eqs. (6.15) and (6.16) to form (6.17), we must 

be aware of the arrangement of the matrices [K.J , j = 1, 2, 3 and [k] t , 

J E-L L 

such that the elements of the matrices can follow the permutation of u., 0., 
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j = 1, 2, 3 in the vector jsj of Eq. (6.17). Furthermore, we notice that, in 
the matrices [p 2 M] and [pB], all elements corresponding to those grids in 
the impact zone are set to zero. Since, as formulated in Section 5, all the 
inertial effects appearing in Eq. (6.7) have been integrated by finite element 
in time a priori, and are included in fK.J 

3 E-L 

It is important to mention that the displacements computed by Eq. (6.17) 
for those grids in the impact zone do not represent the movement of the same 
material particles at any instance. Instead, they represent those for differ- 
ent material particles at various times that occupy the given grids at that 
instance. In the E 'zone, which is composed of material failure, the grid 
points are fixed in the space. In the E-L zone, the grid points move with 
the relative velocity £7 with respect to the particles occupying those grids 
at that instance. Thus, in the E-L zone, the grids must be updated according 
to the definition (Eq. (3.3)) by letting vary from.S7^_ = v^ in the E zone to 
-£7^. = 0 on the interface of the E-L zone and L c zone. 

The solution of Eq. (6.17) for those grids belonging to the Lagrangian 
zone, on the other hand, satisfy the small deformation theory. Thus, the 
grids in this region do not require updating, except those on the interface 
of the E-L and L c zones. It can be expected, however, that the deviation 
on these interface nodes will be negligible compared to those in the impact 
zone. 

6.3 RELEVANT CRITERIA IN CELFE SYSTEM 

We are now ready to formulate the numerical procedure for the com- 
plete CELFE system. The principal part of the CELFE system has been 
concerned with the analysis of the dynamic behavior of the target during and 
right after the penetration process. As mentioned previously, however, the 
complete CELFE system must have the capability to simulate various cir- 
cumstances when certain criteria are detected during the computations. 


6-13 



The criteria related to the CELFE procedure can be classified into 
two categories: (1) the criteria arising from the impact process; and (2) 
the criteria describing the. degeneracy of the. dynamic behavior with respect 
to time due to damping effects. 

The present code consists of three possible impact cases, namely, the 
'penetration process composed of the principal portion of the code; the re- 
bound of the projectile from the target; and the sliding of the projectile on 
the target surface. It is clear that the penetration process follows the failure 
criterion; and the criterion for rebound is rather straightforward. Thus, in 
what follows, only the failure and sliding process will be presented, 

6.3.1 Numerical Procedure for Predicting Material Failure 

The failure criteria for fibrous composites and isotropic materials 
have- been discussed in Section 3. The numerical procedure for predicting 
material failure using these criteria are different in CELFE as presented 
in the following. 

Orthotropic Materials: For fiber composites, Chamis' criterion states 

that: 


F(o^,S^,Kg) >0 no failure 
F (cr^ , ) = 0 incipient failure 

F (°f ’ S f ' ) < 0 failure 


(6.18) 


where 


F^Sj.Kj) = 1 - 


(a 


£ 1 la\ 
11a/ 
2 


/ CT f 22(3 \ J a £ 337 s 


l 22(3/ 


1337/ 


123S \ + ( H 13S 
’i 23S/ \ S fl3Sy 


f 12S 


- K 


i 223 1337 _ 


a illa 337 
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- K 


11a °l 22f3 


iUapS^na S i22p. 


(6.19) 
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(cf. Section 3.4.2). Here, the stresses are computed from Eqs. (6.9) and 
(3.21). 


The application of Eq. (6.19) for the nonlinear analysis of material is 
suggested as follows (Ref. 28): 


• Using the current load cr^ 1 ^ obtained from Eq. (6.9) and Eq. (3.21), 
compute F of Eq. (6.19). 

® Test the value of F according to Eq. (6.18). If F(or^,S,K) = 0, 

then determine the minimum magnitude of the various stress 
ratios: 


Max 


g l 11a 
pi 1 la 


a i22g 

S i22|3 


337 °l 23S 13S °l 12S 

S l 337 ’ S f23s’ S f 13S ? S il2S ; 


where a, (3, and y equal T or C. 

Assuming that the stress ratio (cr^ ZZT^I 22T^ max ^9 aum ’ ^ en 
in the next and subsequent time steps, let E^ ^ V l 23’ 12 

and 23 h e zero. The procedure is similar for other stress 
ratios or combinations. 


isotropic Materials: For isotropic materials, a tentative stress . 
approach suggested in Ref. 10 is used. The procedure can be stated as 
follows: 


* 


© Using the current tentative stresses, (sj^ 1 *^) computed from 
Eq. (6.9), test the yielding criterion. If 


. (n+ 1 ) Z 1 .^(n+l) 


) 



( 6 . 20 ) 


then, the tentative stresses are corrected according to the 
formula 


S. ( ? +1) = (1 - P) (S^ +1) ) 
U ' v i] ' 


( 6 . 21 ) 
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where 


P 


1 - 


yIz/3 Y 2 

[(S-[ +1) ) > (sf M) * 


Wz 


6.3.2 Sliding Due to Oblique Impact 


There are two types of sliding between projectile and target surfaces, 
i.e., the sliding of a projectile on the target surface due to oblique impact, 
and the separations due to spallation. Since the failure zone has been modeled 
by viscoplastic flow, only the sliding of the first type needs to be considered in 
this code. 


Accordingly, the sliding due to oblique impact can be characterized by 
the following two factors: 

• Surface friction between the target and projectile when there is 
no failure on the surfaces; and 

• Slip flow of viscoplastic fluids when failure occurs on the surfaces. 

When there is no failure immediately after the impact, only the first factor 
will be taken into account throughout the sliding process, and the projectile 
will eventually separate from the target. Otherwise, a mixed effect will arise 
and the projectile may either stop and adhere somewhere on the target, or 
rebound away depending on the process. 

The criterion that covers both factors can be obtained using a modified 
Coulomb' s friction law: 


< ll f for no sliding 

(6.23) 

|~tl = jU, l~ni ^° r s ^^ 1Q S onset 

whereof t is the tangential force on the impact surface, ^ is the normal com- 
pressive force on the surface, and jl is a friction factor. For the first case 
where no failure occurs, Eq. (6.19) reduces to the classical Coulomb's law, 



6- 16 



and ju. is the Coulomb friction coefficient. For the second case, can be 

correlated to the slip flow criterion. In the present case, it depends on the 

ratios C 0 . 0 ./C 0 ^._ 0 , i = 1, 2, of the target as well as the projectile, together 
djoj 3 ado 

with the surface tension on the sliding surface where material failure is de- 
tected. In general, we have ji = 1+7, where 7 is the specific free energy 
per unit area. 


6.3.3 Degeneracy Due to Damping 

For large t after the penetration process is completed, the damping 
effects of the structure may smear out the disturbances and reduce the hydro- 
dynamic pressure to some extent, such that the nodal values of the velocity 

v = °C£o) en tire structure (6.24) 

where v is the impact velocity. It can be easily verified that the governing 

'-“O 

equations (6.1) for the impact zone are a system of degenerate equations for 
t > tj, where t^ is some positive constant that Eq. (6.24) holds for all t > t^. 
For < t < t g , where the local acceleration 9v/ 3t remains 0(1), the con- 
tinuity equation and energy equation in (6.1) become trivial; and the momentum 
equations degenerate to (6.3) with' r elation (6.4) after some manipulation. Thus, 
for all t > t^, the dynamic behavior in the entire structure satisfies (6.3), i.e., 
the small deformation theory holds for the entire structure. As t increases 

further to t > t , the damping effects will lead the structure to a static state, 
s 

i.e., Eq. (6.3) will finally approach the Navier equations. 


Based on the formulation given in Section 5, where the conservation 
form Eq. (6.1) was used, the finite element analog will become singular for 
all t > t^. In the present code, a criterion is set to specify the character- 
istic time t^, so that the mesh of the entire system becomes Lagrangian for 
all t > t^. Then, the computations will switch from coupled Euler ian- 
Lagrangian modes to purely Lagrangian satisfying Eq. (6.16) for the entire 
structure. The criterion is set in such a way that the velocity (momentum) 
field and/or pressure field are comparatively uniformly distributed in the 
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entire system (target and projectile) with maximum amplitudes reaching a 
certain fraction of the initial peaks-, and there is no further failure. 

6.4 NUMERICAL PROCEDURE FOR THE COMPLETE CELFE SYSTEM 

In this subsection, the complete procedure for the CELFE system cover- 
ing the previously discussed aspects is presented. The procedure includes 
both the in-core option and the coupled CELFE/NASTRAN option. 


The CELFE procedure can be stated as follows: 

1. The entire CELFE system is integrated with respect to time using 
the generalized two-step time integration scheme discussed in 
Section 6.3. For the coupled CELFE/NASTRAN option, the scheme 
may not be necessary to include the NASTRAN part. In the present 
code, it is included for the sake of simplicity in constructing the 
interface-, procedure. 

2. A unique procedure composed of the following steps is used for both 
the in-core and coupled CELFE/NASTRAN options. For each time 
step: 

• Predict the solutions of primary variables in various zones 
from the corresponding equations using the results obtained 
in the previous time step. 

• Construct the global system matrix equation for the entire 

. structure in terms of coupling variables, and solve the equa- 
tion using the solution obtained in the previous time step. 

• Update the secondary variables in various zones using the 
above predicted values. 

• Update the mesh in the E-L zone and/or L c zone. 

• Correct the results by repeating the process using the 
above predicting solutions together with the results ob- - 
tained in the previous time step. 

• Test the corrected results to observe which criterion 
stated in the Introduction of this section will be satisfied, 
e.g., whether the projectile will rebound, slide, or penetrate, 
etc. 

3. If the test shows the penetration process is still taking place, re- 
peat the procedure stated in (2) for the next time step. Otherwise, 
we have: 

a After the projectile rebounds from the target, and/or 
disturbances in the structure are damped to a certain 


6-18 



fraction of the original values, the program will regard 
the entire structure as a Lagrangian zone and the dynamic 
response will be handled accordingly. 

• If the projectile slides over the target due to oblique im- 
pact, the mesh for the entire structure will be- set to be 
Lagrangian. Then repeat the procedure described in (2). 

The flow chart depicted in Fig. 6-2 illustrates the above stated procedure. 

In order to maintain the accuracy of computations in the vicinity of the impact 
point, as well as a workable in-core storage and a straightforward programming, 
the dimension of the impact zone in this code will be appropriately chosen a prori. 
At t = 0, the entire impact zones are identical with the E-L zone, where the zone 
is initially discretized into Lagrangian mesh, i.e., letting ^ = 0 for all k= 1,2, 3. 
For t > 0, when the failure is detected at certain nodes, the mesh containing these 
failure nodes is turned to Eulerian enveloped by a Lagrangian moving front. Hence, 
the impact zone, whose dimension has to be specified a priori, must be sufficiently 
large to contain the failure region up to the impact process completion. If ex- 
cessively large impact zones are assigned, on the other hand, it will require un- 
necessary computer, storage and computational time. 

6.5 INTERFACING PROCEDURE FOR CELFE/NASTRAN 

The interfacing procedure is constructed for solving large structural 
systems to remove any limitation of the size of the CELFE program due to 
excessive CELFE program core requirements. The procedure is built-in 
as a part of the coupling procedure to handle either or both of the following; 


• Utilizing the matrix equation solver SOLVE of NASTRAN to obtain 
the solution of Eq. (6.9) in the impact zone; and/ or 

® Coupling the Eulerian and Lagrangian modes for the entire struc- 
tural system. 


6.5.1 Utilization of NASTRAN Subroutine SOLVE as Equation Solver 


Symbolically, Eq. (6.9) can be rewritten for each element in the impact 


zone as 


A (a) $ (a) = B (a) 


(6.25) ' 
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Fig. 6-2 - CEEFE Procedure 
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■where a = 1,2,..., M, and M is the total number of the primary variables. 

In this study, M = 14 in general, i.e,, V., j = 1, 2, 3; E; p; and S^., i, j = 1, 2, 3. 

J J 

For orthotropic materials, it reduces to 11, since S-. = S... 

J 


Equation (6.25) can be solved using SOLVE for each a at a time. In 
order to optimize the time required for the interfacing procedure, the follow- 
ing scheme is used: 


/ a \ ( a \ 

® At each time step, compute all A v \ and B'' ' one at- a time for 
each element, and for each a = 1, . . . , M, sequentially.. For 
each a, these element matrices are stored subsequently in 
NAS TRAN INPUTT 2 matrix format on a disk file. 

Store all relevant data, including the solutions of the previous 
time step, and those of the secondary variables, etc., into a 
disk file before calling EXIT from CELFE. 

NASTRAN reads all the data of element matrices A^ and B^, 
and sums the element matrices for each a. Then solve Eq. (6.25) 
using SOLVE for each. 

The solutions a = 1, . . . , M, are written on a disk file using 

OUTPUT2 matrix format. 

• CELFE reads the output disk and makes further computations. 

6.5.2 Coupling Eulerian-Lagrangian Modes Using CELFE/NASTRAN 

In this portion, the main task is to solve Eq. (6.17) for the entire struc- 
ture. Recall that Eq. (6.17) is composed of Eq. (6.15) for the elements in the 
impact zone, and Eq. (6.16) for the elements in the Lagrangian zone. Further- 
more, the element matrices K., G. of Eq. (6.15) and p^M, pB and K of Eq. (6.l6) 

3 J 

for the L c zone are generated and/or computed by CELFE. These element 
matrices are added according to the fixed form (Eq. ( 6 . 1 6 ) ) for the NASTRAN 
part. 








The scheme used in this portion is similar to the previous one with some 
additional steps to incorporate the CELFE substructure into the NASTRAN sub- 
structure as discussed above. The scheme can be stated as follows: 


• Compute all K., j = 1, 2, 3 in Eq. (6. 15) each at a time for each 
3 

element. The element matrices K. for j = 1,2, 3, are stored 

3 
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subsequently on j = 1, 2, 3, for' each element on a- disk file in 
INPUTT 2 format. 

• Compute all G., j = 1,2, 3, and reassemble in a sorted manner 

for all grids of CELFE substructure. Then write G = {G^G^ G^}, 

together with all other data and parameters, in bulk data format 
and punch into card images. These card images together with the 
CELFE grid data are stored in a breakpointed disk file. 

• Store all relevant data in a disk file before calling EXIT from 
CELFE as in Section 6.5.1. 

• NAS TRAN reads all EC, j = 1, 2, 3 from the' disk file and adds the 

matrices into NASTRAN structural stiffness matrix K in Eq. ( 6 . 1 6 ) 
to form the coupled system matrix 

© Add the inertia and damping terms, together with the loading vector 
composed of those read from the breakpointed disk file for CELFE 
grids and the remaining, to form Eq. (6.17). Then solve the system 
and store the results in OUTPUT2 format on a disk file. 

© Return to CELFE for further computations as in Section 6.5.1. 

6.5.3 Remarks 

As mentioned previously, the procedures utilizing SOLVE in NASTRAN 
as an equation solver for the impact zone, and coupling the CELFE substructure 
with NASTRAN can be applied separately or simultaneously. The choice de- 
pends on the specific problem under consideration. The flow chart depicted 
in Fig. 6-3 illustrates a combined version of the procedure. The logic for 
embedding the procedure into the complete CELFE system can be found in 
Fig. 6-2. 
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7. TEST PROBLEMS AND NUMERICAL RESULTS 


In this section, some results obtained from numerical experiments of 
the developed code are discussed. Various test problems have been used in 
various stages of developing the CELFE code. The test problems include an 
inviscid hydrodynamic model, a hydroelasto-viscoplastic model and a complete 
CELFE model. The numerical experiments are carried out to serve the pur- 
pose of developing computer codes stage-by-stage on one hand, and on the other 
hand, to resolve the difficulties in various aspects encountered in the high veloc- 
ity impact analysis. 

7.1 INVISCID HYDRODYNAMIC MODEL IN EULER1AN MODE 

The present problem is used mainly to develop a computer code based 
on the algorithm formulated in Section 5 for handling the viscoplastic flow in 
the impact zone. 

7.1.1 Description of a Typical Impact Problem 

As pointed out in our earlier discussions, a hydrodynamic model is a 
good approximation in the early stages of the high velocity impact process, 
during which the pressure is comparable to the shear strengths of the target 
material. As a starting point we begin with a numerical solution to the hydro- 
dynamic equations with the inviscid adiabatic approximation. In this case the 
conservation equations, with cr.. = -P6.., become 

af ^ = “ &cT ( pv 0 

l 

■k (pv i> = - tr < p > - tr (pv i v j> + 

i 3 

= -ilr< Pv i> - sr- (pev i> +p£ i v i 

i i 


i w (pe > 
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The above equations can be written in the conservation form as 


where 


H 9F k 
S't dx 

k 
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0 = 


P 

V-; j = 1.2,3 
E 


F k = 


v k P 

v k V. + P; j 

v k E+v k P 


= 1,2,3‘ 


; k = 1,2, 3 


(7.2) 
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0 

F ; j = 1, 2, 3 

J 


v.„ F 


i 


; l = l, 2, 3 


(7.3) 


where V. = P v y j = 1, 2, 3, E = pe. 

Consider the impact of a cube on a semi-infinite body composed of pos- 
sibly different materials at a velocity V Q normal to the target surface. A 
typical finite element mesh symmetric about the z-axis is shown in Fig. 7-1. 


Initial Conditions: We assume that the target is at rest at time t = 0. 

Let Sp denote the set of nodal points in the projectile and denote the set of 
nodal points in the target. Then S q = S O denotes the set of nodal points 
common to the target and projectile (i.e., on the interface). For t-= 0, we 
have the following initial conditions for density, pressure, and internal energy: 


Pi 
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(7.4) 
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Fig. 7-1 - A Typical Finite Element Mesh of the Projectile-Target 
Configuration 
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where the subscript "i” refers to the i node, p , and p 

1 ot op 

densities of the target and the projectile, respectively. 


being the initial 


Initial conditions on the velocity components u^, w and w^ in x, y, and 
z-directions, respectively, are 


u. = 
1 

< 

H- 

11 

o 

' ics t ns 

w. = 

o, 

iCS 

X 


P 

w. = 

- V . 

i CS 


(7.5) 


Boundary Conditions: As mentioned earlier the pressure must be zero 
on the free surfaces. Since the target is semi- infinite, the material particles 
far away from the impact region are unaffected. Hence, the internal energy 
and velocities must be zero and the density must be undisturbed. These bound- 
ary conditions, together with conditions for the plane of symmetry, are sketched 
in Fig. 7-2. 


Remark on the' Equation of State; As mentioned earlier, there are two 
forms of equation of state frequently being used in high velocity impact prob- 
lems. They are the Los Alamos equation of state and Tillotson 1 s equation 
of state, both obtained from experiment. Obviously, values of computed pres- 
sure distribution in materials depend heavily on the particular equation of state 
being used, which in turn will affect the entire numerical solution. Some un- 
certainty apparently exists regarding the two equations as they do not appear 
to agree closely in general and each seems to have its own range of validity. 
For instance, the Los Alamos equation of state and Tillotson' s equation of 

state differ substantially in the case of aluminum (see Figs. 7-3 and 7-4) for 

3 

internal energy values below 7.875 Mpa-in /ibm, and agree closely in the 

2 

range of 42 to 45.15 Mpa-in /lbm, as seen in Fig. 7-5. Therefore, caution 
must be exercised when a specific equation of state is employed in the num- 
erical computations. It seems that a more accurate equation of state is 
needed. 
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F.ig. 7-3 - Variation of Pressure with Internal Energy (Aluminum 
P o = 0.0976 lb/in 3 , pc - 0 to 15.0 Mpa) 



Fig. -7-4 - Variation of Pressure with Internal Energy (Aluminum 
- P q = 0.0976 lb/ in pe = 27 to 7,5 Mpa) 
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7.1.2 One-Dimensional Impact Problem 

A one -dimensional impact problem is used herein to test the validity of 
the formulations presented in Section 5. The problem under consideration is 
that an infinite plate traveling at a specified speed, v , hits normally on another 
infinite plate at rest, as depicted in Fig. 7-6a. The two plates are 11.8 in thick, 
and composed of the same aluminum material with initial density P Q = 0,10 lb/in 
The impact velocity is chosen to be v q = 262.5 ft/ sec for which a solution in 
elastic range of impact is available for comparison (Ref. ID). 

The problem was solved by the three-dimensional models described in 

Section 5. The mesh and nodal numbering are generated by subroutine MESH 

and depicted in Fig. 7-6b. The Rankine-Hugoniot pressure, P , is readily ob- 

s 

tained using the Los Alamos equation of state (cf. Appendix), 

P g = 5.99868 x 10~ 2 Mpa (7.6) 

The following numerical computations are performed based on inviscid assump- 
tion together with the Los Alamos equation of state. 

As was discussed in detail in Ref. 37, it is found from the numerical 
solution of Eq. (7.2) that the method of weighted residuals does not perform 
satisfactorily. The pressure development in the material computed by these 
schemes tends to grow indefinitely. This phenomenon becomes more severe 
when the impact velocity is increased. Spurious oscillations with fairly large 
amplitudes always gather near the wave front, and the peak pressure is seen 
to exceed the Hugoniot pressure by about 50' to 100% as the velocity, v, in- 
creases from 262.5 to 24610.0 ft/sec. Figure 7-7 shows the time history of 
pressure at interface computed by the method of least squares. The results 
computed by a modified Galerkin approach using mid-point Runge-Kutta scheme 
are found to be almost identical to those shown in Fig. 7-7. It is seen that the 
pressure at the interface tends to diverge as time increases. The numerical 
instabilities can also be detected by the spurious oscillations near the shock 
fronts as shown in Fig. 7-8. 
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Fig. 7-6 - One -Dimensional Impact Problem 
a. Configuration of Plates at Time 
of Impact 













Pressure, P x 10 {Mpa} 




Fig. 7-8 - Pressure Distributions at Various Times for 1-D Inqpact Problem Using 
Least Squares Scheme with 16 Linear Elements (v = 0.00315 in/psec, 
P q = 0.10 lb/in 3 ) ° 


Solution Computed by Weak Solution Formulation: The pressure history 
at the interface computed using Eq. (5.15) and Eq. (5.16) with a = l/2 is depicted 
in Fig. 7-9, and the corresponding pressure distributions at various times are 
shown in Fig. 7-10. It is seen that the results are improving compared to the 
ones computed by the method of weighted residuals as depicted in Figs. 7-7 and 
7-8. The solution may be improved even further when the parameter a in 
Eq. (5.15) is increased. Figures 7-11 and 7-12 show the pressure history at 
the interface using a = 2 and 4, respectively. The corresponding pressure 
distributions are plotted in Figs. 7-13 and 7-14. 

All results discussed so far and depicted in Figs. 7-7 through 7-14 are 
computed using 16 linear elements. As shown in Figs. 7-10, 7-13 and 7-14, 
the overshoot due to the numerical instability decreases rapidly, and the zig- 
zag behavior near the wave front becomes less severe as the parameter a 
increases. These facts indicate that the numerical dissipation introduced in 
the scheme is proportional to the parameter a. 

The plots also show some substantial differences in amplitudes of the 
pressure waves among different a and a values. In addition, as shown in 
Figs. 7-11 and 7-12, the numerical dissipation increases as a decreases, and 
the phase lag becomes more apparent at the same time. These imply that 
various errors may enter in the computations if the scheme is over -dissipated. 
The sensitivity of the solutions to the factor a, however, can be removed by 
refining the mesh as shown in Figs. 7-15 through 7-18 where 30 even spaced 
linear elements were used .in computations. In this case, the same time step 
size as employed in the l6-element case was used. - This verifies the statement 
mentioned in Section 5.3, namely, the factor a is indeed related directly to the 
ratio of the time step size and the element size. As the mesh size is reduced, 
the effect of a on the accuracy of the solution becomes less important. The 
results using 16. quadratic elements also show the same trend. In Figs. 7-19 
and 7-20, the pressure history at the interface is depicted using a = 2.0 and 
a = 4.0, respectively, with different a. As expected, these results are im- 
proved compared to the cases with 16 linear elements, but are less smooth 
when compared with that using 30 linear elements. Similar comparisons can 
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Fig. 7-10 - Pressure Distributions at Various Times 
with 16 Linear Elements (a = 0.5, a = 0.0) 


30 

'ime, t (jUsec) 


40 


50 


ivelopment at the Interface 
:ar Elements (a = 2.0) 


7-16 



a = 0.00 


a = 0.50 


11 \- 



Fig, 7-12 - Pressure Development at the Interface 
with 16 Linear Elements (a = 4.0) 


7-17 



Pressure, Px 10 (Mpa) 



Fig. 7-13 - Pressure Distributions at Various Times 
with 16 Linear Elements (a = 2.0) 
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Fig. 7-14 - Pressure Distributions at Various Times 
with 16 Linear Elements (a = 4.0) 
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Fig. 7-15 - Pressure Development .at the Interface Using 30 Linear Elements 
{a = 2.0, a = O.'O and 0.25) 
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Fig. 7-16 -Pre ssure Development at the Interface Using 
30 Linear Elements (a = 4.0, a. = 0.0) 


7-21 


Pressure, P x 10 (Mpa) 



•Fig, 7-17 - Pressure Distributions at Various Times Using 30 Linear Elements (a = 2.0, a = 0.0) 
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Fig. 7-18 —Pressure Distributions at Various Times Using 30 Linear Elements (a = 4.0, a = 0.0) 




Fig. 7-19 — Pressure Development at the Interface 
, . Using 16 Quadratic Elements 

{a = 2.0, a = 0.0 and 0.25) 
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Fig. 7-20 — Pressure Development at the Interface 
Using 16 Quadratic Elements 
{a = 4.0. a = 0.0 and 0.25) 
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also be made for the pressure distributions (see, e.g., Figs. 7-13 and 7-14 
for 16 linear elements. Figs. 7-17 and 7-18 for 30 linear elements, and Figs. 
7-21 and 7-22 for 1 6 . quadratic- element's). 


As is seen, all the cases seem to underpredict the peak pressure com- 
pared to the analytic solution. However, with the same given conditions and 

Los Alamos equation of state, the pressure computed by Rankine-Hugoniot 

- 2 

relations is P = 5.99868 x 10 Mpa. This value is approximately equal 
s 

to the average value as shown in Figs. 7-17 and 7-18 which indicates that our 
numerical results are quite reasonable. Accordingly, P g should be the upper 
bound for the average pressure distribution in the material under high velocity 
impact, since the physical as well as numerical dissipative effects may actually 
reduce the pressure buildup in the material. 

» 'j* * 

The speed of sound in a material can be obtained by the relation 



s 


( 9P) 

9P s 


Following the Gibb 1 s equation, one has 


2 _ P 9P 3P 
s " p Z de + dp 


(7.7) 


(7.8) 


where e is the specific internal energy which is now being considered as a 
function of total specific energy and the particle velocity. Thus the propaga- 
tion speed of the pressure waves can be computed by substituting the nodal 
solution of the conservation equations and equation of state on the element 
containing the wave front into Eq. (7.8). The numerical values of Cg in the 
target obtained using Los Alamos equation of state at some typical times are 
shown in Table 7-1. It is seen that the values oscillate around the shock 
speed computed by Rankine-Hugoniot relation, i.e., C~_ TT = 17575.1 ft/sec 
In elas to -dynamic theory, on the other hand, the sound speed is represented by 

C = -J - ■ ^ ~ ^ (7 Q) 

^s Ei >(1 + v) (l-Zv)p U ’^ } 
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Table 7-1 

PROPAGATION VELOCITY OF PRESSURE WAVES IN TARGET 
AT VARIOUS TIMES AFTER IMPACT 
(WITH 30 LINEAR ELEMENTS) 

7 

Poisson Ratio, v = 0.33; Young's Modulus, E = 10 psi; 
Impact Velocity, V q = - 262.50 ft/sec 


Time, t(jU,sec) 

0 

5.0 

10.0 

18.0 

26.0 

31.6 

35.0 

Shock Speed, 
C g (ft/ sec) 

17533.4 

17756.9 

17666.0 

17677.2 

17611.2 

17651.2 

17645.0 


Rankine-Hugoniot — C = 17575. 1 ft/sec 

s 


3-D Elastic Waves — C g = J {1 + ( " x ^ = 19924,2 ft/sec 















■where v is the Poisson ratio and E is Young's modulus. With v = 0.33 and 
E = 10 7 psi for aluminum, Eq. (7.9) gives C s .g^ = 19924.2 ft/ sec, which 
deviates from the values obtained by the present approach by about 10%. 

i 

This implies that the compressibility effects play an important role in the 
dynamic response of materials under high velocity impact loads. In addi- 
tion, although in the present hydrodynamic model the material is assumed 
to be inviseid, the dissipation resulting from the equation of state and num- 
erical viscosity as well may affect the local speed of the pressure waves. 

Momentum and Energy Distributions; Momentum and total energy dis- 
tributions at t = 30 jlsec are depicted in Fig. 7-23 for 16 linear elements. Fig. 
7-24 for 30 linear elements, and Fig. 7-25 for 16 quadratic elements. In these 
plots, the parameters a = 2.0 and a = 0 are used. A better accuracy using 30 
linear elements against 16 quadratic elements is obvious. In particular, as 
seen in the rarefaction region, the momentum and total energy distributions 
computed using l6 quadratic elements are severely distorted by the numer- 
ical instabilities. Hence, a larger parameter a is needed, although it is not 
so obvious from the plots of pressure distributions. 

7.2 HYDROELASTO -VISCOPLASTIC MODEL 

The crucial difficulties encountered in developing the CELFE code are 
caused by the following facts: (1) the occurrence of shockwaves, and/or the 
local existence of large gradients in physical quantities; (2) the coupling 
procedure for global analysis; and (3) the massive computer storage and 
computation time required due to the three-dimensional nature. By ex- 
cluding the latter two factors, the simple one -dimensional inviseid test prob- 
lem discussed in Section 7.1.2 verified the success of the algorithm formulated 
in Section 5 in treating the shock propagation in the impact zone, A 
similar approach is applied in this subsection to demonstrate the actual 
simulation of the impact zone using the weak solution algorithm developed in 
Section 5. 
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Fig. 7-24 - Momentum and Total Energy Distributions at t = 30.0 ^sec 
(30 Linear Elements; a = 2.0, a = 0.0) 
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7.2.1 One -Dimensional Problem in Eulerian Mode 

The one-dimensional impact problem discussed earlier is used also to 
test the present model. The problem description and the finite element mesh 
have been shown in Figs. 7-6a and 7-6b. In the computations, the material 
constants are assumed to be: 

3 

Initial density, p =0.100 lb/ in 
o ^ 

Shear modulus, p = 4.002 x 10 psi 

The impact velocity, v q = 262.5 ft/sec, and Los Alamos equation of state 
are also used for the present problem. 

In order to compare with the test problem in Section 7.1.2, the yielding 
prediction is not incorporated at first. Similar conclusions as discussed in 
Section 7.1.2 can be drawn from the present test case. Figure 7-26 illustrates 
the axial stress and pressure, developments at the impact surface using 30 
linear elements, with the values of a = 4.0 and a = 0.0. 


When the yield strength, Y = 100.5 ksi, is taken into account, as de- 
picted in Fig. 7-27, the stress and pressure buildups decrease due to the 
formation of viscoplastic flow. As time increases, however, the stress and 
pressure approaches the same values' as shown in Fig. 7-26. 

7.2.2 Two-Dimensional Test Case for Unidirectional Fiber Composites — 
Eulerian Mode 

This test case was performed to examine the workability of the code 
for anisotropic materials. The configuration of the problem is depicted in 
Fig. 7-28, with target made of graphite HM/ERLA4617, and projectile of 
aluminum. The material properties are listed in Table 7-2, In these runs, 
the fibers are assumed to be parallel to the x-axis, and the projectile is 
isotropic. 
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Pig. 7 26 - Axial Stress and Pressure Developments at the Interface Using 
30 Linear Elements (a = 4.0, a = 0.0) - Purely Elastic Case 
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Numbering for a Two-Dimensional 
ional Elements) 


Table 7-2 


PHYSICAL PROPERTIES FOR A GRAPHITE 
HM/ERLA46I7 TARGET HIT BY AN ALUMINUM 

PROJECTILE 



Target 

Projectile 

Material 

Graphite 

HM / ERLA4 6 1 7 ' ' . 

Aluminum 

Density 

(lb/in3) 

0.056 

0.101 

Elastic 

Modulus 

(psi) 

E £ll = 27.5 x 10 6 
E i22 = E i33 = 1 '° * 106 

10.5 x 10 6 

Shear- - 

Modulus 

(psi) 

G I12' G «3 = 0 -’ X1 ° 6 
G £Z3 = 0.4 X 10 6 

4.002 x 10 6 

Poisson 1 s 
Ratio 

■ V av 

0.33 

Uniaxial 

Failure 

Stresses 

(ksi) 

s iht = 122,0 • 

S ilic = 128,0 
S ^22T S £33T = 6,1 
S BZC, S £33C = 28,5 
S j?12s** ^‘S. 13s- = 8 * 9 
S £23s = 5,6 

Static Yield: 

40.0 

Fracture Stress: 
.100.5 


7-38 



















Figures 7-29 and 7-30 show the normal stresses and pressure devel- 
oped versus time at nodes 155 and 177, respectively, with an impact velocity, 
v q = 0.0192 cm/jLisec (630 ft/sec), normal to the target surface. The stress 
and pressure distributions at various sections for a given time are illustrated 
in Fig. 7-3 la and b. The onset of the failure in the target is found around 0.06 
jl sec. The failure zone propagates outward uniformly as time increases. 

Runs were also made for v Q = 24610.0 ft/ sec. It was found that, with 
the mesh depicted in Fig. 7-28, the spurious oscillations are severe as shown 
in Fig. 7-32. This indicates that the mesh is too coarse for' an impact velocity 
of this order of magnitude. Failue occurred almost immediately after the 
impact (t = 0.02 jitsec), and the propagation of the failure zone is depicted in 
Fig. 7-33. As shown in Fig. 7-33a, the failure detected at nodes 77, 111 through 
115 may be caused by the numerical instabilities. Thus, finer mesh'is required 
for larger impact velocity. 

7.2.3 Test Case for CELFE Code 

During the development of the CELFE program, various runs of 1-D 
problems were repeated for debugging purposes. The mesh composed of 30 
linear elements with a prescribed impact zone is shown in Fig. 7-34. The 
impact zone contains grids from 29 to 124. Notice that, in this particular 
case, we assigned all the grids to be Eulerian at t = 0, except those interface 
nodes 29 through 32, instead of using Lagrangian mesh (i.e., moving mesh 
with relative velocity £2^ = 0) as suggested in Section 6. Notice also in the 
Lagrangian zone, i.e., the grids 1 through 28, that similar 3-D linear ele- 
ments are used, and can thus be regarded as the L c zone. 

Following the procedure discussed in Section 6, the weak solution 
formulation is used to analyze the impact zone, and the classical approach 
using 3-D solid elements is used to handle the L c zone. The runs were 
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Fig. 7-29 - Stresses and Pressure Developments at Node No. 155 {Target — 

Graphite HM/ERLA4617, Projectile — Aluminum, V q = 630 ft/sec) 
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Fig. 7-30 - Stresses and Pressure Developments at Node No. 177 (Target — 

Graphite HM/ERLA4617, Projectile — Aluminum, V = 360 ft/sec) 
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Fig. 7-31 - Axial Stress and Pressure Distributions at t = 1.0 jusec; (a) on (Y, Z) -Plane at Y = 0, 
Z = 3.0315 in.; (b) on'the Axis of Symmeti f (i.e., X = Y = 0) (Target — Graphite HM/ 
ERLA46 17, Projectile — Aluminum, V Q = 630 ft/sec) 
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Fig. 7-32 - Axial Stress Distribution at Time t = 0,2 ^ sec 
(V = 24610.0 ft/ sec) 
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Fig. 7-34 - Finite Element Mesh of CELFE for 1-D Problem 
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made using the in-core CELFE option. The numerical results are illustrated 
by the axial stress distributions at various times in Figs. 7-35 and 7-36. The 
distributions are compared with the results solved by the Eulerian code as dis- 
cussed in Sections 7.1 and 7.2, and the Lagrangian mode constructed using the 
classical approach for the zone., 

The results show that the CELFE solution is almost identical to the 
Eulerian solution at all times. 

7.3- CEEFE ANALYSIS OF THREE-DIMENSIONAL IMPACT PROBLEM 

The results of a 3-D impact problem are discussed in this subsection to 
demonstrate the CELFE in-core as well as the CELFE/NASTRAN operation. 
The given condition of the example can be summarized as follows: - 

Target: Boron AVCO 5505 plate with dimensions: 

9.75 in. x 8.75 in. x 0.0’82 in. 

. ~ ‘ 3 3 

Projectile: Silastic with dimensions: (0.30) in 

Impact 

Velocity: 630 ft/sec normal to the plate 

The geometric configuration is depicted in Fig. 7-37, and the physical prop- 
erties for the materials are listed in Table 7-3. Tbe boron epoxy. target is 
assumed to be unidirectional with fibers parallel to the y-axis. 

7.3.1 Numerical Results for CELFE In-Core Runs 

Due to the symmetry, only a quadron of the plate is cons.idered. The 
finite element mesh of the model is shown in Fig. 7-38, where the impact 
zone consists of the first 33 elements (nodes 1 through 76), with elements 
26 through 33 (nodes 62 through 76) as the translation zone (E-L zone). The 
remaining are assigned to the L c zone. 

The computations show that the failure starts at 0.001 /tsec after the 
impact. The failure front propagates as the time increases as illustrated 
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Fig. 7-35 - Axial Stress Distributions at t=5.0fisec Using 30 Linear Elements 
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Fig. 7-36 - Axial Stress Distributions at t = 12.0 fj . sec Using 30 Linear Element's 


Fig. 7-37 


Sketch of Three-Dij 
Unidirectional Bore 
Silastic Particle 


msional Impact Problem of a 
AVCO 5505 Composite by a 





Table 7-3 

' PHYSICAL' PROPER TIES FOR A BORON AVCO 5505 
PLATE HIT BY A SILASTIC PROJECTILE 



Target 

Projectile 

Material 

Boron AVCO 5505 

Silastic 

Density 

(lb/in3) 

0.073 

0.0512 

Elastic Modulus 
(psi) 

E m = «.2xio 6 
E i22 = 33 = 3,15 x 10 

• 

Shear Modulus 
(psi) 

■ G I1'2 •= G il3 =' 0-78-xlO 6 

G«3 = °- 6xl ° 

0.6 x 10 b 

Poisson's Ratio 

"az = "ns = °-”' 

v tZ 3 = °' 53 


Uniaxial Failure 
Stresses (ksi) 

3 illT = ^99*‘0 

s mc ‘ 

1 

S i22T = S i33T = 8 ' 1 

S i22C , = S i33C “ - 17 ' 9 

S *12S = S il3S = 9,1 

S i23S = 8 - 9 

Yield; 

2.7027 
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Fig. 7-38 - CELFE Finite Element 
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in Fig. 7-39. The pressure development versus time and the normal stress 
distributions are shown in Figs. 7-40 and 7-41, respectively. 

7.3.2 CE LFE/NAS TRAN Runs 

The NASTRAN program adopted in the present development has been 
the 15.7 version from the NASA Lewis Center. During the development, the 
erroneous INPUTT2 and OUTPUT2 modules were corrected, and some sub- 
routines were modified to make the interfacing procedures more efficient. 

The detailed discussions will be presented in Part II of this report. 

In order to compare the results computed by CELFE in-core runs and 
CELFE/NASTRAN runs, the size of the elements for both cases are made 
similar. In the CELFE/NASTRAN model, the finite element mesh for the 
impact zone is identical to the CELFE in-core case, and is depicted in 
Fig. 7-42. It consists of the first 33 linear isoparametric elements. In 
the Lagrangrangian substructure, the L c zone is assigned to be empty for 
simplicity, and all elements in the substructure are composed of'NASTRAN 
general quadralateral elements, i.e., L^ zone. The mesh is illustrated in 
Fig. 7-43. To ensure the displacement compatibility across the boundary, 
multipoint constraints (MPC) are employed to connect the CELFE and ! 
NASTRAN substructures. The connecting procedure can be illustrated as 
follows: A NASTRAN grid point is added to each CELFE mid-point grid 
lying on the E-L zone, which has now become the interface of the CELFE 
and NASTRAN substructure due to the empty L^ zone, e.g., the grid's 163, 

166, 169, 172 and 175 in the NASTRAN substructure (c.f.. Fig. 7-43) are - 
added, respectively, to the CELFE grids 63, 66, 69, 72 and 75 (Fig. 7-44). 

The mid-point grids are coupled by MPCs for the three translational degrees 
of freedom; rotational displacements are transmitted from the quadralateral 
plate elements to the CELFE linear isoparametric elements by using rigid 
bars connected between the mid-point node, say, node 163 depicted in Fig. 
7-44, and the top and bottom nodes of the solid, nodes 64 and 62, respectively. 

The numerical results computed by the CELFE/NASTRAN model are 
found to be nearly identical to the CELFE in-core model. The deviation 
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Fig. 7-39 - Propagation of Failure Front with Respect to Time 
(Nodes shown are at mid-plane.) 
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Fig. 7-40 - Pressure Developments at Various Nodal Points (Fig. 7-39) 
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Fig. 7-42 - CELFE Finite E! 
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Fig. 7-43 - NASTRAN Finite Element Model in CEEFE/NASTRAN Runs 
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between the runs is of 0 (10 ). The spurious oscillations found in the 

CELFE/NASTRAN results, however, are slightly greater than the in-core 
runs. 

7.3.3 Concluding Remarks 

As mentioned previously, the numerical results obtained by CELFE 

in-core and CELFE/NASTRAN runs are almost identical, although in the 

fomal case, 3-D linear isoparametric elements are used in the Lagrangian 

zone, and the later, 2-D plate elements are used. For the computing time, . 

on the other hand, it requires approximately six minutes, per time step for 

the in-core operation, in contrast to 30 minutes for the CELFE/NASTRAN 
% 

operation. The substantial increase of run-time in CELFE/NASTRAN opera- 
tion arises mainly due to the interfacing procedures between the CELFE and 
NASTRAN programs; yet the generality of the NASTRAN program may as well 
add some inefficiency in the computations. 

During .the CELFE/NASTRAN executions, it is found that the high-core 
is always required in the NASTRAN transient. analysis module. For a 10- 
element (with '7 elements in the impact zone) model, 70K of core was required 
in the Univac 1108 system; for 16 elements, about 90K of core was set. In 

t » 

the present model of 54 elements (with 33 elements in the impact zone), 100K 
of core was assigned to the NASTRAN program. Apparently, there are in- 
appropriate overlays and/or other errors in the transient module of the 15.7 
version. This problem also added some inefficiency in the CELFE/NASTRAN 
executions. 

i 

It is noted that, in the present demonstration problem, the thickness 

of the plate is much less than the width of the plate; the ratio of the thickness 

-2 

to width (in both x- and y-axis) is e = 0.082/5, = 0 (10 ). In the impact zone, 

the finite element mesh consists of elements with a thickness of 0.041 in., and 
width (in both x- and y-axis) varies from 0.15 in. to 1.0 in., i.e., the ratio varies 
from = 0 (10 *) to = 0 (10 ^). Although isoparametric elements are used 
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in the impact zone where the mesh is mapped to unit cubes represented in 
the local coordinate system. With 1 j the finite element approximation 

of the gradient of a physical quantity, say, <f>, with respect to the global co- 
ordinate system, may still induce an error that may be magnified by 
0(l/^F h ). As a consequence, the mesh with ratio « 1 may lead to the 
possibility of imposing a singular problem in the finite element approxima- 
tion. This phenomenon is precisely the boundary layer type of problem en- 
countered in-dealing with high Reynolds number flows. 

This phenomenon was demonstrated by the present problem in various 
numerical tests using various mesh arrangements. According to the experi- 
ment results, conducted by I. M. Daniel and T. Diber [38], there is no failure 
under the impact conditions specified in this test problem. However, the 
numerical results computed using various mesh arrangements show that the 
failure occurs almost immediately after the impact. It is also shown that, 
for the 10-element model, (where the average mesh ratio = 0.041/0.15) 
the failure front propagates faster than the other models with more regular 
mesh (including the in-core models), even though some measures were made 
to smear out the spurious oscillations. The results are improved when more 
elements in the impact zone become more regular, i.e., with ratio = 0(1). 

The discrepancy of the numerical results from the experiment data 
may partially be induced also by some other factors, such as the accuracy 
of the equation of state and the failure criteria. It is clear that these factors 
are equally or even more important than that of the singular mesh used in the 
impact zone to cause the inaccuracy in the computations. Hence, further 
investigations on this aspect is necessary in order to improve the accuracy 
of the CELFE program. 
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8. SUMMARY AND DISCUSSION 


In the present study, a three-dimensional, coupled Euler ian-Lagrangian 
finite element (CELFE) program was developed for high velocity impact. The 
study was divided into two phases. Phase I was devoted to constructing a 
numerical algorithm for analyzing the local dynamic behavior in the vicinity of 
the impact point during the penetration process. The local analysis was based 
on a hydroelasto-viscoplastic model governed by conservation equations of 
mass, momentum and energy, an appropriate set of constitutive equations, 
and an equation of state. A computer code based on the theorem of weak 
solutions was developed in terms of the Eulerian mode. 

In Phase II, the CELFE program was developed to extend the results 
obtained in Phase I to perform the global analysis of certain structures sub- 
jected to high velocity impact. The computer program was composed of 
three parts: the first part was generalized from the code developed in Phase I 
for handling the local analysis in the impact zone represented by a general 
moving coordinate system instead of the Eulerian system; the second- part was a 
classical Lagrangian mode for handling the dynamic response of the remaining 
part of the structure induced by the impact; the third part was a coupling pro- 
cedure to integrate the above two parts. 

The coupling procedure, in turn, consisted of two parts. The first part 
was a procedure to smoothly translate the coordinate system, and thus the 
dynamic representation, from the Eulerian mode describing the dynamic be- 
havior in the failure zone to the Lagrangian mode describing the dynamic re- 
sponse away from the impact zone. The procedure was incorporated in the 
CELFE program by the variation of the relative velocity with respect to the 
local velocity of material particles. The second part of the coupling pro- 
cedure was constructed for interfacing the CELFE and NASTRAN programs 
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to remove any limitation of the size of the CELFE code due to excessive 
CELFE core storage requirements. The interfacing procedure consisted of 
several I/O routines installed in CEEFE to read- and write the matrices in 
NAS TRAN INPUTT2 OUTPUT 2 formats. The large matrix equations for the 
conservation and constitutive equations, together with the data for the CELFE 
substructure to be coupled with the NASTRAN substructure, were then carried 
into the NASTRAN program to perform all the required manipulations using 
the standard NASTRAN Rigid Format 9. 

During the development, various test cases were performed to investi- 
gate the validity and capability of the algorithm in various aspects. The nu- 
merical experiments show that the most difficult part in analyzing a structure 
subjected to high velocity impact, namely, the highly nonlinear, large deforma- 
tion, and shock and failure propagation phenomena in the vicinity of the impact 
point, can be resolved using the finite element algorithm developed in Section 5. 
The massive computer core storage required for three-dimensional analysis 
has also been partly resolved by coupling the Eulerian and Lagrangian modes, 
and interfacing with NASTRAN (or any other large structural program). 

As this is only the first stage of the development, the CELFE program 
will need further improvements and extensions. One of the immediate tasks 
would be a further reduction of the in-core storage by improving the arrays 
of existing variables in the program. 

From the numerical experiments, it is found that the results will be 
sensitive to the equation of state. In addition, the effects of possible phase 
changes due to impact can be included via an appropriate equation of state. 

Thus, further studies on the thermodynamic aspects of the high velocity im- 
pact problem are necessary. Along this path, additional thermal effects can 
be incorporated into the program by modifying the governing equations. 
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Finally, as stated in Section 7.3.3, a substantial increase of run time 
(compared to the in-core execution) is needed in performing the CELFE/ 
NASTRAN computations. This fact shows the inefficiency of using the inter- 
facing procedure. Hence, for the long run, it is suggested that the CELFE 
program should be incorporated in the NASTRAN as a module. 
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RANKINE-HUGONIOT RELATIONS 



Appendix ' 


During the short period of time in which the projectile impacts on the 
target, a plane shock wave is generated in the projectile as well as in the 
target. Due to the extremely high pressure behind the shock front, the ma- 
terials can essentially be assumed to behave as inviscid, compressible fluids. 
The problem, including the geometric development of the impact shock, can 
then be treated as an unsteady, supersonic flow resembling a moving shock. 


Assuming that the hemispherical shockwave profile is steady in time, 

the Rankine-Hugoniot relations relating the pressure, P, internal energy, 

and density, p, behind the shock to the same quantities in front of the shock 

are applied. These equations express the conservation of mass, momentum, 

and energy in terms of the shock velocity, v , and particle velocity, v : 

s P 

P v = p (v - v ) {A. 

'o s s p 


P - P 

o 


P V V 
o s p 


(A. 2) 


|(e - e ) - v 2 /2 p v = P v 
o p'J os os 


(A. 3) 


where the subscript o refers to the initial (or undisturbed) values. The pro- 
duct p_ v is called the shock impedance. Using (A.l) and (A. 2), (A. 3) can be 
o s 

written as 


e - e o = i< p + p o 

o 


Solving for v and v from (A.l) and (A. 2), we obtain 
s p 


/ t (p -p) n 

v s= y V P = Vp o> ^o- 1 ) 


(A.4) 


(A. 5) 


A-l 



The initial pressure, P , is generally very small and can be neglected in 
comparison with P. If shock velocity and particle velocity are known at a 
point on the' shock front, the pressure .can be -computed- -us mg {A. 2). One 
more relation is needed to solve the flow across a shock, namely, the equa 
tion of state; a specific form of the equation of state is given in Subsection 
3.4. 
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